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(D . Abstract 

Q 

This is the second in a series of papers describing a 3+1 computational scheme 
^N ■ for the numerical simulation of dynamic black hole spacetimes. In this paper 

pq ' we focus on the problem of numerically time-evolving a given black-hole- 

CN . containing initial data slice in spherical symmetry. We avoid singularities 

vQ \ via the "black-hole exclusion" or "horizon boundary condition" technique, 

^D \ where the slices meet the black hole's singularity, but on each slice a spatial 

neighbourhood of the singularity is excluded from the domain of the numerical 
computations. 
Q^. We first discuss some of the key design choices which arise with the black 

Vh . hole exclusion technique: Where should the computational domain's inner 

boundary be placed? Should a free or a constrained evolution scheme be 
used? How should the coordinates be chosen? 
^ ' We then give a detailed description of our numerical evolution scheme, 

j^ ■ We focus on a standard test problem, the time evolution of an asymptotically 

flat spherically symmetric spacetime containing a black hole surrounded by 
a massless scalar field, but our methods should also extend to more general 
black hole spacetimes. We assume that the black hole is already present 
on the initial slice. We use a free evolution, with Eddington-Finkelstein-like 
coordinates and the inner boundary placed at a fixed (areal) coordinate radius 
well inside the horizon. 

Our numerical scheme is based on the method of lines (MOL), where the 
spacetime PDEs are first discretized in space only, yielding a system of coupled 
ODEs for the time evolution of the field variables along the spatial-grid-point 
world lines. These ODEs are then time- integrated by standard finite difference 
methods. In contrast to the more common "space and time together" finite 
differencing schemes, we find that MOL schemes are considerably simpler to 
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implement, and make it much easier to construct stable higher order differ- 
encing schemes. Our MOL scheme uses 4th order finite differencing in space 
and time, with 5 and 6 point spatial molecules and a 4th order Runge-Kutta 
time integrator. The spatial grid is smoothly nonuniform, but not adaptive. 

We present several sample numerical evolutions of black holes accreting 
scalar field shells, showing that this scheme is stable, very accurate, and can 
evolve "forever" . As an example of the typical accuracy of our scheme, for a 
grid resolution of Ar/r ~ 3% near the horizon, the errors in gij {Kij) compo- 
nents at t = 100m in a dynamic evolution are ^ 10~^ (3 x 10~^) in most of the 
grid, while the energy constraint is preserved to ~ 3 x 10~^ of its individual 
partial-derivative terms. 

When the black hole accretes a relatively thin and massive scalar shell, 
for a short time 3 distinct apparent horizons are present, and the apparent 
horizons move at highly superluminal speeds; this has important implications 
for how a black-hole-exclusion computational scheme should handle the inner 
boundary. 

Our other results for the scalar field phenomenology are generally consis- 
tent with past work, except that we see a very different late-time decay of the 
field near the horizon. We suspect this is a numerical artifact. 
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I. INTRODUCTION 

Dynamic black hole spacetimes are interesting both as laboratories in which to study the 
basic physics and phenomenology of strong-field gravitation, and as probable very strong 
astrophysical sources of gravitational radiation. These systems are strongly relativistic, 
time-dependent, and often highly asymmetric, and so are hard to study in detail except by 
numerical methods. 

This is the second in a series of papers describing a 3 + 1 computational scheme for the 
numerical simulation of dynamic black hole spacetimes. In the first paper in this series 
(Thornburg (1999), hereinafter paper I), we described our numerical initial data construc- 
tion, and presented some example initial data slices. In this paper we discuss our numerical 
time-evolution scheme and related topics, and present several example spacetimes computed 
using this evolution scheme. 

Black hole spacetimes necessarily contain singularities, which a 3 + 1 computational 
scheme must somehow avoid in order for the 3 + 1 equations to remain well-defined and 
numerically tractable. Traditionally this avoidance has been done by using "freezing" slic- 
ings, where the lapse function is forced to drop to a very small value in the spatial vicinity 
of a singularity-to-be, causing proper time to cease to advance there and thus "freezing" 
the evolution before it reaches the singularity. Unfortunately, by advancing at very different 
proper-time rates in different parts of space, such slices necessarily "stretch"^, causing the 
3 + 1 field variables to develop increasingly steep gradients and non-smooth features, usually 
near to or just inside the horizon. This then leads to a variety of numerical problems (see, 
for example, Smarr (1984); Shapiro and Teukolsky (1986)). 

To avoid these problems, Unruh (1984) suggested evolving only the region of spacetime 
outside the black holes. This "black hole exclusion", or "black hole excision", or "horizon 
boundary condition" technique^ has been developed by a number of researchers (e.g. Thorn- 
burg (1985, 1987, 1991); Seidel and Suen (1992); Thornburg (1993); Scheel et al. (1995a); 
Anninos et al. (1995b); Scheel et al. (1995b); Marsa (1995); Marsa and Choptuik (1996); 
Scheel et al. (1997); Cook et al. (1998)). This technique involves choosing slices which in- 
tersect the singularity, but on each slice excluding (excising) some spatial neighbourhood 
of the singularity from the domain of the numerical computations. With singularities thus 
avoided, the slicing may then be (and generally is) chosen to avoid or minimize slice stretch- 
ing throughout the computational domain. The boundary of the excluded region, the "inner 
boundary" of the computational domain, is typically placed on or somewhat inside the 
apparent horizon. 



^The term "grid stretching" is often to describe this, but we prefer "slice stretching", as the 
effect is actually inherent in the continuum freezing-slicing 3 + 1 equations, and is unrelated to the 
numerical grid. 

^We prefer the first two terms as focusing attention on the underlying process - the exclusion of 
the black hole from the computational domain - rather than on the particular boundary condition 
used for implementation. Indeed, in our computational scheme the boundary in question is actually 
well inside the horizon, so "horizon boundary condition" would be something of a misnomer. 



Using this technique, a well-designed evolution scheme can evolve an initial slice "forever" 
without encountering coordinate singularities, thus attaining the "Holy Grail of numerical 
relativity" (Shapiro and Teukolsky (1986, page 76)), "... a code that simultaneously avoids 
singularities, handles black holes, maintains high accuracy, and runs forever." 

In this paper we first discuss some general aspects of the black hole exclusion technique, 
focusing on the treatment of the inner boundary and on how the coordinates should be 
chosen. We then give a detailed description of our black-hole-exclusion numerical evolution 
scheme for dynamic black hole spacetimes. Finally, we present several sample spacetimes 
computed using this evolution scheme. 

Our goal is a computational scheme which can be generalized to a wide range of black 
hole spacetimes, so as much as possible we first present and develop our algorithms for the 
generic case of a spacetime with no symmetries and arbitrary matter fields present, and 
only later specialize to a specific testbed system. However, the numerical results discussed 
here all pertain to a relatively simple testbed system: an asymptotically fiat spherically 
symmetric spacetime containing a black hole surrounded by a (massless, minimally coupled) 
scalar field. ^ 

The spherically symmetric scalar field and similar systems have been studied by a number 
of past researchers, including (among others) an extensive analytical study by Christoudolou 
(1986a,b, 1987a,b, 1991, 1993), other valuable analytical work by Gomez and Winicour 
(1992a), 3 + 1 studies by Choptuik (1986, 1991); Seidel and Suen (1992); Choptuik (1993b); 
Bernstein and Bartnik (1995); Scheel et al. (1995a,b); Anninos et al. (1995b); Marsa (1995); 
Marsa and Choptuik (1996), 2 + 2 studies by Goldwirth and Piran (1987); Goldwirth et al. 
(1989); Gomez and Winicour (1992b); Hamade and Stewart (1996), a hybrid 3 + 1 and 2 + 2 
study by Gomez et al. (1996), and an interesting comparison of 3 + 1 and 2 + 2 methods 
by Choptuik et al. (1992).'^ Of these, our work is closest in spirit to that of Marsa (1995); 
Marsa and Choptuik (1996), but we use slightly different spatial coordinates, a very different 
inner-boundary treatment, and very different numerical methods. 

In the remainder of this paper, we first summarize our notation (section II), then dis- 
cuss the tradeoffs of where the inner boundary should be placed in a black-hole-excluded 
evolution scheme (section HI), whether a free or a constrained evolution scheme should be 
used (section IV), and various general criteria for choosing the coordinates (section V). We 
then discuss our formulation of the continuum 3 + 1 geometry and scalar field equations, our 
coordinate conditions, and our boundary conditions for the evolution and coordinate equa- 
tions (section VI). We then discuss our numerical methods (section VII) and our diagnostics 
for assessing the numerical results (section VIII). We then present and analyze a number of 



^We have previously developed and implemented a similar scheme for the evolution of a vacuum 
axisymmetric spacetime containing a single black hole (Thornburg (1993)), but for that system 
we were unable to obtain time evolutions free of finite differencing instabilities. We believe these 
instabilities were primarily due to the handling of the z axis in the polar spherical coordinate 
system, and were not intrinsic to the black hole exclusion technique, but we haven't proven this. 

^There has also been a wide variety of work (both analytical and numerical) by many authors 
investigating the self-similar behavior discovered by Choptuik (1993b). 



sample numerical evolutions (section IX). We end the main body of the paper with some 
conclusions and directions for further research (section X). Finally, in the appendices we 
tabulate the detailed evolution equations for the spherically symmetric scalar field system 
(appendix A), give a brief introduction to the method of lines (appendix B), summarize 
our methodology for convergence testing and briefly review the technique of Richardson 
extrapolation (appendix C), and describe our techniques for trying to estimate the effects 
of floating-point roundoff errors by artificially adding a tiny amount of noise to the field 
variables during an evolution (appendix D). 

II. NOTATION 

We generally follow the sign and notation conventions of Misner et al. (1973), with a 
(—,+,+,+) spacetime metric signature, G = c = 1 units, and all masses and coordinate 
distances also taken as dimensionless. We assume the usual Einstein summation convention 
for all repeated indices, and we use the Penrose abstract-index notation as described by (e.g.) 
Wald (1984). However, for pedagogical convenience we often blur the distinction between a 
tensor as an abstract geometrical object and the vector or matrix of a tensor's coordinate 
components. 

We use the standard 3 + 1 formalism of Arnowitt et al. (1962) (see York (1979, 1983) for 
recent reviews). For our spherically-symmetric-scalar-field test system, we use coordinates 
{t,r,6,(f), with 6 and ip the usual spherical-symmetry coordinates. However, in general we 
leave t and r arbitrary, i.e. unless otherwise noted we make no assumptions about the specific 
choice of the lapse function or the radial component of the shift vector. 

The distinction between 3- and 4-tensors is usually clear from context, but where am- 
biguity might arise we use prefixes ^^^ and ^'^^ respectively, as in ^^^R and ^^^R. Any tensor 
without a prefix is by default a 3-tensor. £^ denotes the Lie derivative operator along the 
4- vector field v"'. 

We use abed for spacetime (4-) indices, and da denotes the spacetime coordinate par- 
tial derivative operator d/dx"'. Qab denotes the spacetime metric and Va the associated 
4-covariant derivative operator. 

We use ijkl for spatial (3-) indices, di denotes the spatial coordinate partial derivative 
operator d/dx"^. gij denotes the 3-metric of a slice, Vj the associated 3-covariant derivative 
operator, and g the determinant of the matrix of Qi/s coordinate components, a and /?* 
denote the 3 + 1 lapse function and shift vector respectively, n"' denotes the (timelike) 
future pointing unit normal to the slices. Kij = ^^nQij = —^i^j denotes the 3-extrinsic 
curvature of a slice, and K = Ki^ its trace, p = n"'n'^Tah and j* = —UaT"-'^ denote the locally 
measured energy and 3-momentum densities respectively. Tij denotes the spatial part of 
stress-energy tensor, and T = Tj* its trace. 

diag[- ■ ■] denotes the diagonal matrix with the specified diagonal elements. x{6)y denotes 
the arithmetic progression x, x + 6, x + 26, x + 36, ■ ■ ■ , y. Gaussian (x= A, a=B) denotes 
the Gaussian function exp(— ^2;^), where 2; = (x — A)/ B. For a symmetric rank 2 covariant 

tensor Tij, we define the pointwise "magnitude" tensor norm {{Tij)) = JTijT"^^. 

When discussing finite difference molecules, we often refer to a molecule as itself being 
a discrete operator, the actual application to a grid function being implicit. We denote 
a molecule's "application point", the grid point where the molecule's result is defined, by 



underlining the corresponding molecule coefficient. We define the "radius" of a molecule in 
a given direction to be the maximum distance in grid points away from the application point 
in that direction, where the molecule still has a nonzero coefficient. We define a "centered" 
molecule as one having the same radius in both directions. For any coordinate x, Ax denotes 
a uniform-grid finite difference grid spacing in x. Thus, for example, we describe the usual 
3 point centered 2nd order 2nd derivative molecule by 
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and say that it has radius 1 in both directions. 

Consider a family of molecules all approximating the same continuum (differential) op- 
erator, one molecule in the family being centered (with radia tq in both directions), and the 
remaining molecules being off-centered. If an off-centered molecule has radia r_ and r_|_ in 
the — and -|- directions respectively, we formalize the intuitive notion of "off-centering dis- 
tance" by defining its "offset" to be (the positive integer) tq — r_ if r_ < tq, or (the negative 
integer) r^ — tq if r^ < tq (we assume that exactly one of these cases holds). Finally, we 
define the offset of a centered molecule to be 0. Thus, for example, we say that the 4 point 
off-centered 2nd order 2nd derivative molecule 
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(Ax) 
which has radia r_ = and r+ = 3, has offset +1. 
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III. THE LOCATION OF THE INNER BOUNDARY 

In 3 + 1 numerical relativity we generally don't know the spacetime we're simulating in 
advance, but rather we construct it "on the fly" during a numerical evolution. In using the 
black hole exclusion technique, then, the precise region of spacetime to be excluded from 
the numerical evolution must in general be chosen dynamically in each slice. ^ In this section 
we consider how this choice - specifled in practice by the inner boundary location - should 
be made. 

We assume throughout this discussion, and indeed everywhere in this work, that a black 
hole or holes is already present in our initial data, and hence that every slice contains a 
black hole. (This restriction could probably be weakened, but we haven't investigated this.) 



^Choptuik (1992a) has suggested that a multipass evolution algorithm might be useful here, 
with the future behavior of the black holes being approximately known from a previous lower- 
resolution evolution of the same spacetime. This would ease the "lack of advance knowledge" 
problem considerably, but an "on the fly" scheme would still be needed for the initial evolution of 
a previously unknown spacetime. 



A. Event or Apparent Horizon? 

The excluded region must clearly satisfy three conditions: 

1. There must be no singularities in the non-excluded "exterior" region of spacetime (the 
computational domain), i.e. any singularities must be contained within the excluded 
region of spacetime. 

2. Furthermore, to prevent numerical difficulties there must be no singularities "close" 
to the exterior region, i.e. any singularities must be at least some (strictly positive) 
distance away from the computational domain, this distance being uniformly bounded 
from below during an evolution. 

3. The evolution equations on the exterior region of spacetime must be well-posed 
(causal), i.e. no future pointing null geodesic may cross from the excluded region 
to the exterior region. In other words, everywhere on the inner boundary the light 
cones must point exclusively out of the computational domain and into the excluded 
region.^ 

In terms of these conditions alone, the ideal choice for the inner boundary location would 
be the event horizon. Condition 3 would then hold by definition, and condition 1 would be 
precisely the well-known cosmic censorship hypothesis, which is already widely assumed 
in 3 + 1 numerical evolutions. Unfortunately, the event horizon can't be located during 
a numerical evolution, as it's defined (Hawking (1973); Hawking and Ellis (1973)) in an 
inherently acausal manner: it can only be determined if the entire future development of the 
slice is known. (As discussed by Anninos et al. (1995a); Libson et al. (1996), in practice the 
event horizon may be located to good accuracy given only the usual numerically generated 
approximate development to a nearly stationary state, but the fundamental acausality still 
remains.) 

Instead, we assume that our spacetime and slicing are such that each slice contains one 
or more apparent horizons, one of which we use as a reference point for placing the inner 
boundary. We first consider the simplest case, where the inner boundary coincides with the 
outermost apparent horizon. (For pedagogical convenience, in this and the next subsection 
we assume there's only a single black hole present, but this restriction isn't in any way 
essential to our arguments.) 

Condition 3 then follows from the definition of an apparent horizon, which, in the words 
of Hawking and Ellis (1973, page 323), ". . . moves outward at least as fast as light; and 
moves out faster than light if any matter or radiation falls through it."^ Israel (1986b, a) has 



^ Using a common analogy with hydrodynamics, we might imagine trying to simulate fluid flow 
in the exterior region, and requiring that fluid should be flowing into the excluded region at a 
supersonic speed. 

^Hawking and Ellis (1973) don't prove this assertion, but it's easily demonstrated (Unruh and 
Wald (1993)) by considering the expansion of the congruence of outgoing null geodesies orthogonal 



also shown that under suitable technical conditions, an apparent horizon on an initial slice 
can always be extended in time in a "rigid" (locally area-preserving) manner so as to act 
". . . as a wall that permanently seals off its interior contents from causal influence on the 
environment." However, in general this extension differs from the time development of the 
apparent horizon. It would be interesting to investigate this extension's use in black-hole- 
exclusion inner boundary placement, but so far as we know this hasn't yet been done. 

With the inner boundary coinciding with the apparent horizon, condition 1 is equivalent 
to requiring all singularities to lie within (outermost) apparent horizons. In terms of generic 
black hole spacetimes, this "apparent cosmic censorship" assumption is quite strong, in 
fact considerably stronger than "standard" cosmic censorship. In particular, Wald and Iyer 
(1991) have shown that even in Schwarzschild spacetime there exist (angularly anisotropic) 
slices which approach arbitrarily close to the r = singularity, yet contain no apparent 
horizons. (The related work of Abrahams et al. (1992b) is also relevant here.) 

In practice, we simply follow the usual course in numerical relativity of assuming all 
necessary properties of spacetime - in our case including apparent cosmic censorship (in 
fact the somewhat stronger condition 2) and the light cones pointing entirely inwards on the 
inner boundary - and on this basis proceeding with the numerical computations, verifying 
our assumptions after the fact. 

B. On or Inside the Apparent Horizon? 

Our analysis of the previous subsection takes the inner boundary to be located coincident 
with an apparent horizon,^ typically the outermost one, and this location has been used 
successfully by Thornburg (1993); Marsa (1995); Marsa and Choptuik (1996); Scheel et al. 
(1997). Alternatively, the inner boundary may be placed somewhat inside an apparent 
horizon (again typically the outermost one), as has been done successfully by Seidel and 
Suen (1992); Scheel et al. (1995a); Anninos et al. (1995b). These latter authors all chose to 
leave a macroscopic distance (i.e. one large compared to the grid spacing) between the inner 
boundary and the nearest apparent horizon, but still another option would be to reduce this 
distance to only a few grid points, i.e. to place the inner boundary a few grid points inside 
an apparent horizon. (We don't know of any numerical tests of this last option.) 

Placing the inner boundary coincident with an apparent horizon has several advantages: 



to an apparent horizon. By definition, this expansion is zero on the apparent horizon. But the 
Raychaudhuri equation (more precisely its analog for null geodesies) implies that the derivative 
of this expansion along the congruence is negative semidefinite, so the congruence must be non- 
expanding on future slices. This congruence must therefore lie within the trapped region in these 
(future) slices, and hence lie on or within their (outermost) apparent horizons. 

*If a multidimensional Cartesian-topology grid is used, then the inner boundary necessarily has 
an irregular "staircase" shape. However, if this is chosen to approximate the apparent horizon to 
within a single grid spacing, then for present purposes we ignore the non-smoothness, and still 
refer to such an inner boundary as "coincident" with the apparent horizon. 



• As discussed in section III A, it ensures the well-posedness (causality) of the exterior 
region evolution (at least so long as each slice does indeed contain an apparent horizon). 
(In contrast, if the inner boundary is somewhat inside an apparent horizon, then the 
causality of the exterior region evolution depends on the details of the spacetime, the 
slicing, and the spatial coordinate choice.) 

• It avoids "wasting" any grid points by placing them inside an apparent horizon, where 
they'd be unable to causally influence the exterior evolution. (In practice, though, this 
inefficiency doesn't seem to be serious.) 

• It's relatively easy to implement, particularly if the coordinates and finite difference 
grid use a polar spherical topology, as then the coordinates and grid may be chosen so 
the relevant apparent horizon is (i.e. coincides with) a coordinate and grid sphere. 

• The condition that the initial data have an apparent horizon coinciding with the 
inner boundary, provides a natural inner boundary condition for the York initial data 
algorithm, and also makes easy to use this algorithm to construct multiple-black-hole 
initial data (Thornburg (1985, 1987)). 

• The condition that the relevant apparent horizon remain coincident with the inner 
boundary during the evolution, provides a natural inner boundary condition for the 
radial component of the shift vector. In contrast, if the inner boundary is inside 
the horizon then enforcing such a condition is more difficult, as the corresponding 
"boundary condition" still needs to be enforced at the horizon, despite this no longer 
being a boundary of the numerical grid. 

However, placing the inner boundary somewhat inside an apparent horizon (whether by 
a macroscopic distance or only by a few grid points), also has a number of advantages: 

• It helps to prevent any numerical inaccuracies at the boundary from affecting the 
exterior evolution. Thornburg (1991); Seidel and Suen (1992); Thornburg (1993); 
Scheel et al. (1995a); Anninos et al. (1995b) all cite this reason for placing the inner 
boundary inside the horizon in their black-hole-excluded computational scheme. 

• A black-hole-excluded evolution often requires that the apparent horizon or horizons 
be explicitly located, i.e. that the apparent horizon equation be numerically solved, 
at each time step. Like the exterior evolution, it's desirable for the apparent-horizon- 
finding process to be unaffected by any computational errors at the inner boundary. 

• Because the apparent horizon equation can't be solved exactly in closed form, numeri- 
cal apparent-horizon finders invariably use iterative algorithms. (See, e.g., Baumgarte 
et al. (1996); Thornburg (1996); Anninos et al. (1998) for recent discussions of nu- 
merical apparent-horizon-finding algorithms.) These algorithms require the 3-1-1 field 
variables {gij and Kij) to be defined throughout a spatial neighbourhood of the appar- 
ent horizon being located. Placing the inner boundary somewhat inside an apparent 
horizon naturally provides this. In contrast, placing the inner boundary coincident 
with an apparent horizon would require some form of extrapolation to continue the 
field variables inside the "inner boundary" for the horizon finder's use. This seems 
inelegant, and might also be numerically ill-conditioned. 

9 



• In a dynamic spacetime where we don't know the apparent horizon's motion in ad- 
vance, the coordinate position of the inner boundary and/or the value of the shift vec- 
tor there must in general be dynamically adjusted to "track" the apparent horizon's 
motion, i.e. to maintain the apparent horizon at a desired (possibly time-dependent) 
coordinate position as the evolution proceeds. As discussed by Thornburg (1993, sec- 
tions 3.1.5-3.1.7 and 4.5), depending on the algorithms used, this "horizon tracking" 
may be somewhat imprecise, i.e. at times the inner boundary may drift somewhat in 
its position relative to the apparent horizon. If this is the case, and the drift happens 
to move the inner boundary outwards, then the numerical evolution may cease to be 
well-posed (i.e. it may violate the causality condition 3 of section III A) unless there's 
a "buffer zone" between the apparent horizon and the inner boundary. 

Considerable further research is needed to investigate these issues, and to help clarify 
the relative merits of the different possible inner boundary placements. 



C. Multiple Apparent Horizons 

The problem of placing the inner boundary is more complicated when a slice may contain 
multiple apparent horizons, either enclosing distinct volumes or topologically nested inside 
one other. In terms of our discussion in section III A, it's generally preferable to place the 
inner boundary near (on or somewhat inside) the outermost apparent horizon within any 
given black hole. This minimizes the chances of any singularities being too close to the 
exterior region (and also minimizes the number of grid points "wasted" by being inside the 
outermost apparent horizon and thus causally isolated from the exterior evolution), while 
still ensuring that the exterior region evolution remains well-posed. 

In practice, multiple apparent horizons are often present only for a limited time interval 
in an evolution. As an example of the general type of phenomenology to be expected, 
consider the 400.pqwl evolution discussed in section IX C, where the initial slice contains a 
black hole surrounded by a relatively thin and massive scalar field shell. Figure 7 shows the 
apparent horizon positions as a function of coordinate time t for this evolution. Notice that 
initially there is only one apparent horizon, then at t ~ 19.13 two new apparent horizons 
form outside the original one, then at t ^ 19.71 the original (inner) and middle apparent 
horizons merge and disappear, leaving only the (new) outer apparent horizon at later times. 

Because the original apparent horizon disappears at t ^ 19.71, a black-hole-exclusion 
evolution scheme which had previously been using it as a reference for placing the inner 
boundary, must somehow switch to using the new outer apparent horizon at later times. 
[Similarly, in a black-hole-collision simulation, when the two black holes coalesce and their 
previous (separate) apparent horizons disappear, a black-hole-exclusion evolution scheme 
must change from using their separate apparent horizons as inner-boundary reference posi- 
tions before the coalescence, to using the resulting single black hole's new apparent horizon 
afterwards (assuming it has one).] 

There are a number of ways such a transition might be made. For example: 

• One possibility would be to smoothly move the inner boundary from the inner ap- 
parent horizon out to the outer one during the time interval (19.13 ^ t ^ 19.71 in 
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our example) where the shces contain both. Unfortunately, this requires a substan- 
tial inner-boundary motion, possibly at a speed faster than the coordinate speed of 
light, or even faster than the grid CFL limit. (For our example, a coordinate distance 
Ar ^ 3.0 must be covered in a coordinate time At ^ 0.6, so some part of the motion 
must necessarily be at a speed in excess of 5 times the coordinate speed of light.) Such 
rapid motion may well cause numerical problems, both near the inner boundary and 
possibly (if the coordinate and grid conditions don't sufficiently attenuate the inner 
boundary's motion) also further out in the grid. At the very least, such highly "super- 
luminal" inner-boundary motion would likely require some type of causal differencing 
(Alcubierre (1993); Alcubierre and Schutz (1994)). 

• To avoid individual grid points having to move so far and so fast, one could instead 
simply drop (delete) from the grid a suitable number of the innermost grid points at 
some suitable time. This would have the effect of discontinuously moving the inner 
boundary outwards, while leaving the spatial coordinate system and the remainder of 
the grid unchanged. 

• A variant of this latter possibility would be to again drop some of the innermost grid 
points, but simultaneously recoordinatize the slice and/or (discontinuously) change 
the grid points' spatial coordinate positions, so as to maintain some desired horizon- 
tracking condition, such as there being a grid point or layer of grid points precisely on 
the horizon. Note that for the black-hole-coalescence case the computational domain 
changes topology as well as size and shape, so some discontinuous change in the grid 
is essential. 

As with the initial placement of the inner boundary, the treatment of such multiple- 
apparent-horizon slices is greatly complicated if spacetime - and in particular the behavior 
of the apparent horizons - isn't known in advance, but rather is computed only "on the 
fly" during a numerical evolution. (For example, in the context of figure 7, however the 
transition from the original to the new outermost apparent horizon is made, it's probably 
desirable to do it somewhat before the original apparent horizon disappears. But how can 
this impending disappearance be detected in advance?) 

Again, much further research is needed to properly investigate these and other possibil- 
ities for how such multiple-apparent-horizon spacetimes and slicings should be treated. 

D. Our Present Scheme 

In our present computational scheme we take a very simplistic approach to these issues: 
we place the inner boundary at a fixed time-independent spatial coordinate position (radius) 
^ = ^miii) chosen to be well inside the horizon's coordinate position at all times. Our code 
explicitly computes the light cones on the inner boundary, and monitors them to ensure the 
exterior region evolution remains causal.^ 



^In detail, the check is that the outgoing light-cone speed c+ defined by (13), satisfies c+ < 
on the inner boundary at all times. This condition has always been satisfied throughout all our 
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This inner boundary placement has proved adequate for our present purposes, but clearly 
wouldn't suffice for more general spacetimes, particularly (non-spherically-symmetric) ones 
containing moving black holes. 

IV. FREE VERSUS CONSTRAINED EVOLUTION 

3 + 1 evolution schemes are customarily classified into "free" , "partially constrained" , 
or "fully constrained" evolution schemes, depending on whether none, some, or all of the 
constraints are used in evolving the field variables. 

It was once thought (e.g. Piran (1983); Stewart (1984)) that the finite differenced con- 
straint and evolution equations are inherently inconsistent, rendering free and constrained 
evolution schemes fundamentally different in meaning, and calling into question whether 
either can truly yield good approximations to solutions of the continuum 3 + 1 equations. 
However, it later became clear from the work of Choptuik (1986, 1991) that there is in fact 
no such inconsistency, i.e. that properly constructed free and constrained evolution schemes 
are both valid, and (if they use the same order of finite differencing) yield results which are 
within an 0(1) factor of each other in accuracy. 

With one (important) exception discussed below, the choice of free versus constrained 
evolution is thus purely a pragmatic one, to be made on grounds of finite differencing sta- 
bility, easy of implementation, and 0(1) factors in accuracy and efficiency. 

There is one serious problem with constrained evolution schemes, though: They generally 
entail elliptic or parabolic (constraint-incorporating) updating equations for the gij and 
Kij components. These equations require boundary conditions, and while outer boundary 
conditions pose no (conceptual) problem, with the black-hole-exclusion technique there's no 
obvious way to obtain inner boundary conditions during an evolution. Further research to 
develop constrained black-hole-exclusion evolution schemes would be very desirable. 

In this work we have chosen a free evolution, both because of the inner-boundary- 
condition problem, and because free evolution schemes are easier to implement and easier 
to generalize efficiently to multiple spatial dimensions. 

V. GENERAL CONDITIONS FOR THE COORDINATES 

How should the coordinates be chosen in a black-hole-exclusion evolution? To permit 
easy interpretation and accurate numerical solution of the 3 + 1 equations, we suggest requir- 
ing the coordinates to satisfy the following general "well-behavedness" criteria throughout 
the computational domain (including near to and on both the event and apparent horizons): 

1. The slices must be spacelike and asymptotically flat. 

2. The coordinates must cover all the "interesting" regions of spacetime (certainly all of 
spacetime outside the black hole) to the future of the initial slice. 



evolutions. 
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3. The coordinates must be nonsingular (smooth) and nondegenerate (spacetime-filhng) 

4. Similarly, the light cones should be non-degenerate, and in particular they should not 
be "infinitely closed up" or "infinitely wide open" as they are (for example) near the 
horizon in the Schwarzschild slicing of Schwarzschild spacetime. 



5. All the 3 + 1 field variables should be smooth functions of the coordinates. 
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To allow well-behaved long-time evolutions, we also require these conditions to continue 
to hold in the limit t — *> oo. (Requiring condition 5 in this limit specifically rules out the 
type of sharp gradients and non-smooth features which arise from slice stretching.) Finally, 
we'd like our coordinate conditions to be readily generalizable to generic multiple-black-hole 
spacetimes. 

Our conditions for the coordinates are very similar to those of Shapiro and Teukolsky 
(1986, page 76), who describe the "Holy Grail of numerical relativity" as "... a code that 
simultaneously avoids singularities, handles black holes, maintains high accuracy, and runs 
forever." 

Schwarzschild spacetime provides a particularly useful test case for considering coor- 
dinate choice. As a number of previous researchers have observed (e.g. Marsa (1995); 
Marsa and Choptuik (1996)), (ingoing) Eddington-Finkelstein coordinates {t,r,6,(f)^^ 
are ideally suited for use with the black-hole-exclusion technique. Figure 1 shows the 
Eddington-Finkelstein slicing of Schwarzschild spacetime, plotted in both Kruskal-Szekeres 
and Eddington-Finkelstein coordinates. These coordinates satisfy all of our criteria 1-5 de- 
scribed above. Note that the Eddington-Finkelstein slices are not maximal: K is nonzero, 
and in fact spatially variable, throughout the slices. 

For example, the light cones shown in figure 1(b) are well-behaved near to and on the hori- 
zon, whereas a similar plot in Schwarzschild coordinates (Misner et al. (1973, figure 32.1(a))) 
shows pathological behavior there. More generally, in these coordinates the 4-metric (and 
hence all the 3-1-1 field variables) remains well-behaved everywhere in spacetime - including 
near to and on the horizon - except near the singularity itself. 

In terms of our discussion in sections III B and III D, it's instructive to check the causality 
of the exterior region evolution for the combination of Eddington-Finkelstein coordinates 



^^This implicitly assumes that the underlying spacetime geometry (e.g. the 4-Riemann tensor) 
is itself smooth. Note that if the type of self-similar behavior discovered by Choptuik (1993b) is 
present, this requirement means we must exclude at least a neighbourhood of any self-similarity 
points. 

^^Recall (Misner et al. (1973, box 31.2)) that these coordinates are defined by taking r to be 
an areal radial coordinate, 6 and (p the standard angular coordinates, and choosing the slicing so 
that t + r is an (ingoing) null coordinate. For reference, we list some of the 3 + 1 field variables 
for Schwarzschild spacetime in these coordinates in appendix A of paper I; Thornburg (1993, 
appendix 2) gives a more extensive list. 

^^Recall that x{6)y denotes the arithmetic progression x, x + 6, x + 26, x + 36, ■ ■ ■ , y. 
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(slices) and an r = constant inner boundary position. For this case it's easy to see that 
the exterior region evolution does in fact remain causal for any (stationary) inner boundary 
position inside the horizon. This can be seen as the slopes of the outgoing legs of the light 
cones plotted in figure 1(b): these point inwards for r < 2m and outwards for r > 2m. 

For Kerr spacetime, Kerr coordinates^^ are similarly well suited for use with the black- 
hole-exclusion technique, and similarly satisfy all of our criteria 1-5 described above. The 
Kerr-coordinates 4-metric is similarly well-behaved everywhere outside the ring singularity, 
including near to and on the (event) horizon. 

For Kerr slicing and a stationary r = constant inner boundary position, the exterior 
region evolution remains causal provided the inner boundary remains between the inner and 
outer horizons (Thornburg (1993, section 3.1.5)). This restriction isn't a problem so long as 
the black hole has non-maximal spin, which we require in any case for the validity of our 
computational scheme. [An maximal-spin (extremal) Kerr spacetime is infinitesimally close 
to a spacetime with no event horizon and a naked singularity, rendering our fundamental 
black-hole-exclusion technique invalid.] 

VI. CONTINUUM PHYSICS 

We now discuss the basic continuum physics of our testbed system and our computational 
scheme. 



A. 3 + 1 Geometry and Scalar Field Equations 

1. Generic Spacetime 

We first briefly review our formulation of the continuum 3 + 1 equations in their generic 
form, making no assumptions about spacetime's symmetries. 

As is customary for the 3 + 1 formalism, we assume that spacetime is globally hyperbolic, 
and we introduce global spacetime coordinates x" = (t,x*), where t is a global time coor- 
dinate. Given the spacetime metric ^^^gab, "we then define the usual 3-metric gij, (3-) lapse 
function a, and (3-) shift vector jS"^, so the spacetime line element becomes 

^^Us^ = ^^^Qab dx" dx^ = -{a^ - A/3') dt^ + 2A dt dx' + g^j dx' dx^ . (3) 

Given gij, we define the usual 3-covariant derivative operator Vj and (3-) Ricci tensor and 
scalar Rij and R = Ri\ 

Taking n" to be the future pointing timelike unit normal to the slices, we define the usual 
(3-) extrinsic curvature Kij = \£ngij = —^i^j. Given the spacetime stress-energy tensor 



^^Recall (Misner et al. (1973, box 33.2)) that these generalize Eddington-Finkelstein coordinates 
for Schwarzschild spacetime: the Kerr slicing in Kerr spacetime is defined by taking the usual 
(areal r), and requiring that t + r be an (ingoing) null coordinate. Thornburg (1993, appendix 2) 
lists some of the 3 + 1 field variables for Kerr spacetime in these coordinates. 
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Tab = ^^^Tab, we define the usual 3 + 1 scalar field variables Tij = ^^^Tij, p = n""n''Tab, and 
f = ~naT''\ 

The 3 + 1 evolution and constraint equations are then (York (1979, 1983)) 

dtQij = - 2aKij + £^gij (4a) 



dtKij = - V.Vj-a + a [Rij - 2g''K,kKji + KKij + £pK, 



^v 



+ A-Ka ((T - p)g,, - 2Ti,) (4b) 

£pgij = V^ + VjA (4c) 

= P'^dkgij + gjkdjP'' + gjkdjp'' (4d) 

£pKij = P'^VkKjj + K,kVj(3' + KjkV^(3'' (4e) 

= /^'^Sfeir,, + K.kd.f^'^ + K,fca,/5'= (4f) 



and 



C=[R- KijK'^ + K' j - (I67rpj = (5a) 

C' = iVjK'^ - V'K) - (Sttj*) = (5b) 



respectively, where the shift vector Lie derivative terms are underlined for reasons discussed 
in section VII B. We use a free evolution scheme, i.e. after constructing initial data (discussed 
in detail in paper I) we use the constraints only as diagnostics of our computational scheme's 
accuracy. 

We take the scalar field to satisfy the 4-scalar wave equation VaV"0 = 0, and to have 
the usual stress-energy tensor 

47rT„b = {da(t>){db(t>) - \gab{dc(t>W(t>) ■ (6) 

We then define the 3 + 1 scalar field variables 

(7a) 
(7b) 

u. ^ 
so that 

(8a) 
V,0 = Pi . (8b) 

These latter definitions (7) and (8) are similar, but not identical, to those of Choptuik (1986, 
1991, 1993b); Marsa (1995); Marsa and Choptuik (1996). 

Straightforward but tedious calculations then give the scalar field evolution equations as 



P^ 


= V.0 




Q 




- /9*V,0 




dt<l) = aQ + (3'Pk 



^tP^ = V,[aQ + P'Pk) (9a) 

dtQ = Vi [aP') + aKQ + P'ViQ (9b) 
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and the 3 + 1 scalar field variables as 

Anp = \ (PkP" + g' 



Anj, = -P,Q 
ATiTi, = P,P, + \g,j (-PkP" + Q' 



47rT = -^PkP'' + IQ^ 



(10a) 
(10b) 
(10c) 
(lOd) 



For the black hole exclusion technique, we assume that each slice contains one or more 
apparent horizons. As discussed by (e.g.) York (1989), an apparent horizon satisfies the 
equation 



H = ViTi' + Ki.n'n^ -K = 0, 



(11) 



where we take n* to be the outward-pointing unit normal to the horizon, and for future use 
we define the "horizon function" H as the left hand side of this equation. 

2. Spherical Symmetry 

We now assume that spacetime is spherically symmetric, and that the spatial coordinates 
are x* = (r, 6, ip), with the usual polar spherical topology. We take the 3 + 1 field tensors to 
have the coordinate components 

(12a) 
(12b) 
(12c) 
(12d) 

Notice that we do not factor out either a conformal factor or any r^ factors from the 3-metric 
components. ^^ 

With these assumptions it's then straightforward, though tedious, to express all the 
other 3 + 1 variables in terms of the "state variables" A, B, X, Y, P, and Q, and their 
spatial (1st and 2nd) derivatives. For reference, we list the resulting equations for the 
constraints and most of our diagnostics in appendix B of paper I; we list the equations 
for our remaining diagnostics in section VIII A of this paper, and those for the evolution 
equations in appendix A. 

Note that although spherical symmetry means there can be only one black hole in space- 
time, there may be multiple apparent horizons. We use h to denote the or an apparent 



9ij = ( 


liag 


A B Bsin^e' 


Kij = diag 


X Y YsinH' 


13'^ 


' f3 O" 




P^ = 


P 








^^ Although this somewhat simplifies the 3 + 1 equations, it may degrade the accuracy of our 
computational scheme. In particular, it's at least plausible that factoring out r^ factors might give 
improved accuracy at large r. It might be interesting to try a side-by-side accuracy comparison 
between otherwise-identical factored and unfactored schemes, but we haven't investigated this. 
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horizon's coordinate position (radius). We describe our methods for numerically locating 
apparent horizons in section Bl of paper I. 

Given our functional form for the metric, it's easy to see that the range of causal spatial 
velocities within the light cone is given by 



B. Coordinate Conditions 

Although we have experimented with a number of different slicings and spatial coordinate 
conditions (some of these are described in Thornburg (1993)), in this paper we discuss only 
a single, relatively simple, slicing and spatial coordinate condition: 

As suggested by Choptuik (1993a); Marsa (1995); Marsa and Choptuik (1996), we choose 
the lapse function a by following the Eddington-Finkelstein requirement that t + r be an 
(ingoing) null coordinate. This gives the "generalized Eddington-Finkelstein" slicing condi- 
tion 

" +/3 = 1. (14) 



'A 

We choose the shift vector /9* to keep the areas of constant-r surfaces (coordinate spheres) 
temporally constant during the evolution, i.e. dtQee = 0, giving the "constant area" spatial 
coordinate condition 

- 2aY + pdrB = . (15) 

This shift vector condition generalizes the usual areal radial coordinate: if the radial coordi- 
nate on the initial slice is areal {geg = r^), then the constant-area coordinate condition will 
preserve this throughout the evolution. 

[In practice we normally do choose our initial data to have an areal radial coordinate, 
but this is solely for convenience in data analysis; our computational scheme doesn't require 
it. Indeed, in our computational scheme the time evolution commutes with 3-coordinate- 
transformations, so a different spatial coordinate choice on the initial slice would leave the 
slicing and the spatial-coordinate world lines invariant, merely relabelling the latter. Such 
invariance is a useful property of a computational scheme (Thornburg (1993, section 2.6)).] 

These coordinate conditions are easy to implement and very cheap to to compute, re- 
quiring only a few algebraic operations at each grid point to compute a and j3. Moreover, 
they satisfy all of our general coordinate-choice criteria described in section V. Unfortu- 
nately, these coordinates implicitly require the existence of a "radial" coordinate, so they 
don't generalize well to multiple-black-hole spacetimes. Moreover, the constant area spatial 
coordinate condition can only constraint a single degree of freedom of the shift vector at 
each event, which is insufficient to fully specify the coordinates in a nonspherical spacetime. 
We will discuss more general and generalizable coordinate conditions in a future paper. 

In anticipation of future generalization to different coordinate conditions, we specifically 
do not assume any particular coordinate conditions in deriving any of our equations outside 
this section, i.e. elsewhere in this paper we take the lapse function and shift vector as 
generic. That is, for example, despite the constant-area shift vector condition (15), we still 
numerically evolve B = ggg in the same manner as the other state variables. 
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C. Boundary Conditions 

Conceptually, our computational domain should be all of spacetime exterior to the inner 
boundary r = rmin- However, as is customary in 3 + 1 numerical evolutions of asymptotically 
flat spacetimes, we further truncate our computational domain at a finite outer boundary 
radius r = r-y. 



max- 



1. Evolution Equations 

As discussed in section III, provided 0+ < on the inner boundary, the evolution equa- 
tions need no (continuum) boundary conditions there. (We discuss the finite difference 
treatment of the inner boundary in section VII B.) 

We impose approximate outgoing radiation boundary conditions at the outer bound- 
ary. This is most elegantly done by matching the numerically computed field variables to a 
tensor-spherical-harmonic decomposition of the gravitational field, as discussed by Thorne 
(1980); Abrahams and Evans (1988); Abrahams (1989); Abrahams and Evans (1990); Abra- 
hams et al. (1995).^^ However, such boundary conditions are quite complicated to derive and 
implement. Instead, here we use a much simpler approximation, the simple pseudo-scalar 
Sommerfeld-type outgoing radiation condition of Trautman (1958a,b); Fock (1959).^^ This 
condition has been used successfully in 3-dimensional (x, y, z) numerical relativity calcula- 
tions by Nakamura and Oohara (1989). 

For our formulation of outgoing radiation boundary conditions (an adaptation of that of 
Thornburg (1993)), we let u denote any of the state-vector field variables A, B, X, Y, P, or 
Q, and u^g denote the same field variable at the corresponding coordinate position in some 
"background" slice to be chosen later. We first factor out u and Wbg's asymptotic coordinate 
dependence at large radia by defining a constant parameter m such that that u = 0{r"^) at 
large radia in a fiat spacetime. 

We then define 6u = u — -Ubg and take 5u/r^ to (approximately) satisfy a simple flat- 
spacetime scalar outgoing radiation boundary condition 

5u fir — c+t) , , _^, 

-^^ ^ - (at r = r,nax) (16) 



0-! 



lyi 1 1 L lyi t I 



for some (unspecified) function /, where c+ is the outgoing light-cone speed defined by (13), 
and n is a constant parameter chosen based on the expected asymptotic fall-off of outward- 
propagating perturbations in 6u. 



^^See also the gravitational-radiation-theory survey of Pirani (1962). 

^^Givoli (1991) gives an excellent survey of numerical methods for non-reflecting boundary condi- 
tions in various areas of computational physics. As well, Israeli and Orszag (1981) have suggested 
a "sponge filter" technique which may be useful in improving the accuracy of such outer boundary 
conditions; Choptuik (1986); Marsa (1995); Marsa and Choptuik (1996) have used this technique 
successfully in numerical relativity computations similar to ours. 



Taking the time derivative of (16), assuming the background field variables to be time- 
independent, and simplifying, we obtain the approximate scalar outgoing radiation boundary 
condition 



df 6u ~ — c_i 



dr Su H Su 



(at r = rmax) • (17) 



All our results in this paper use the outer-boundary-condition parameters given in ta- 
ble I, and take the background slice to be an Eddington-Finkelstein slice of a mass-mtotai 
Schwarzschild spacetime, where mtotai is the total Misner-Sharp mass of the current slice. ^^ 
We discuss the efficacy of these boundary conditions in section IX F. 

2. Coordinate Conditions 

Our coordinate conditions discussed in section VI B are algebraic equations in a and 
f3, and hence need no boundary conditions. However, for consistency with more general 
coordinate conditions which we will discuss in a future paper, we use outer (though not 
inner) boundary conditions for our coordinate conditions here. 

As discussed by York (1980), we assume a ^ 1 and /?* ^ at spatial infinity to ensure 
that our coordinates are asymptotically Minkowskian. 

For the lapse function, we approximate this at our finite outer boundary radius with the 
scalar Robin outer boundary condition suggested by York and Piran (1982), 

(V,a)n* + ^^ = O (^^) ^ (at r = r^^x) • (18) 

York and Piran (1982) describe a vector Robin condition which properly takes into 
account the shift vector's tensor character, but for simplicity, at present we instead just use 
a scalar Robin condition independently for each coordinate component of the shift vector, 

(V./5^>^ + ^ = O (^-^) ^ (at r = r^ax) • (19) 

In spherical symmetry, these conditions become 

-^ + ^- = O (- j ^ (at r = r_) (20) 

and 



^^^ + - = O[i^]^0 (at r = r^aj (21) 



A r \r 



respectively. 



^^mtotah ^^^d hence the background field variables, are actually slightly time-dependent, violating 
one of our assumptions in deriving (17). This effect is generally small, and we ignore it. 

19 



VII. NUMERICAL METHODS 

We now describe our methods for discretizing and numerically solving the 3 + 1 evolution 
and coordinate equations. 

The fundamental assumption underlying all our numerical methods is that all the 3 + 1 
field variables are smooth functions of the spacetime coordinates throughout the compu- 
tational domain. That is, our numerical methods are specifically not designed to handle 
shocks, gravitational chaos, self-similar behavior of the type found by Choptuik (1993b), 
or any other phenomenology which generates significant power at arbitrarily high spatial 
frequencies. (This smoothness requirement is the motivation of our coordinate-choice crite- 
rion 5 in section V.) 

A. Nonuniform Gridding 

To allow adequate resolution of the rapidly changing field variables near the black hole 
while still extending the grid well out into the weak-field region, some sort of nonuniform 
grid spacing is highly desirable. 

An adaptive gridding scheme such as that of Gropp (1980); Berger (1982); Berger and 
Oliger (1984); Berger and Jameson (1985); Berger (1986) is probably the most efficient way 
to do this, and Choptuik (1986, 1989, 1992b, 1993b); Masso et al. (1995); Briigmann (1996) 
report excellent results in using such schemes in spherically symmetric 3 + 1 black hole 
evolutions. However, such adaptive gridding schemes are quite complicated to implement. 

For present purposes we have instead chosen a much simpler non-adaptive nonuniform 
gridding scheme, where the grid spacing varies smoothly in a predetermined (and in our case 
time-independent) manner across the computational domain. For the (fairly limited) class 
of spacetimes we deal with here, such a scheme provides a substantial efficiency gain over a 
simple uniform grid, but is still easy to implement. 

To describe our nonuniform gridding scheme, we introduce a new radial coordinate w 
such that the grid is uniform in w. Although it would be possible to use (w, 6', ip) coordinates 
as a tensor basis as well as for the numerical grid (see, example, Abrahams et al. (1992a); 
Bernstein (1993)), we prefer to retain the natural {r,6,ip) coordinates for this purpose, 
treating the nonuniform gridding as a lower-level numerical technique transparent to our 
continuum problem formulation. 

We therefore implement nonuniform gridding by rewriting individual r-coordinate partial 
derivatives in terms of w-coordinate partial derivatives, i.e. we write 

dr = J(9^ (22a) 

a,, = J^a^^ + Kd^ , (22b) 

where we define the coefficients J = dw/dr and K = d'^w/dr'^. 

It remains to specify how the nonuniform grid coordinate w should be chosen. The log- 
arithmic coordinate w = log(r/ro) is often used (e.g. Abrahams et al. (1992a); Bernstein 
(1993)), but we find that this offers insufficient control over how the grid resolution varies 
as a function of (spatial) position. In particular, with such a logarithmic coordinate we find 
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that a grid giving good resolution near the black hole, will typically have insufficient reso- 
lution to accurately represent asymptotically-constant-width features in the field variables 
(e.g. outward-propagating scalar field shells) at large radia. 

To obtain finer control over how the grid resolution varies as a function of position, we 
instead use the "niixed-210" nonuniform grid coordinate (Thornburg (1993, section 7.1.7)) 
defined by 

dw 1 11 

dr a(r/ro)^ 6(r/ro) c ' 

where Tq, a, b, and c are suitable (constant) parameters. (The name "mixed-210" comes 
from dw/dr being a mixture of 1/r^, 1/r^, and 1/r^ terms.) Integrating and taking w = 
at r = ro, we have 



- = - 1-^ +xlog - +- --1 • (24) 




We also need the inverse function r{w); we compute this numerically.^^ 

The physical interpretation of this coordinate is that for suitable choices of a, b, and c, 
in the innermost part of the grid, the first term in (23) is dominant, so the grid spacing 
is roughly Ar = a(r/ro)^ ■ Aw. In the middle part of the grid, the second term in (23) is 
dominant, so the grid spacing is roughly Ar = 6(r/ro) ■ Aw. In the outermost part of the 
grid, the third term in (23) is dominant, so the grid spacing is roughly Ar = c ■ Aw. 

By convention, we always place the inner grid boundary at w = 0, i.e. r = ro = rmin. 
All our results in this paper use the parameters ro = 1.5, a = oo (effectively omitting the 
first term in (23) and (24)), b = 5, and c = 100. Figure 2 shows w{r) and Aw{r) for these 
parameters; w qualitatively resembles a logarithmic radial coordinate in the inner part of the 
grid, and a uniform radial coordinate in the outer part of the grid. In particular, note that 
the grid spacing Ar asymptotes to a constant value at large radia, allowing asymptotically- 
constant-width features to be accurately evolved even in the outer parts of the grid. 



B. Spatial Finite Differencing 

Our spatial finite differencing uses a single non-staggered grid, uniform in w, with grid 
points located at the fixed spatial coordinate positions w = w^i^{Aw)wraa.x- We label the 
grid points with the integer grid coordinates w = Winin(l)wmax- 

To construct a finite difference approximation to a continuum equation, we first apply 
the nonuniform-gridding transformations (22) to replace all r-coordinate partial derivatives 



^^We use Brent's ZEROIN subroutine (Forsythe et al. (1977); Kahaner et al. (1989)) to numer- 
ically find a zero of the function w{r) — w^,, where ti;* is the (given) value of w for which r is 
desired. This technique is easy to program and very robust and accurate, but mildly inefficient as 
it doesn't exploit the smoothness of w{r). However, in practice we only need r{w) at grid points, 
so we simply precompute and store all the r values. 
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with tf-coordinate partial derivatives. (Or equivalently, we could use (22) to numerically 
compute the r-coordinate partial derivatives in terms of the w-coordinate ones.) We then 
discretely approximate each w-coordinate partial derivative with a suitable finite difference 
molecule. 



1. Grid Interior 

In the grid interior, we use the usual centered 5 point 4th order molecules for most (9^ and 
dujyu terms in the 3 + 1 equations. However, for the (9^ terms derived from the shift vector 
Lie derivative terms in the 3 + 1 evolution equations (these terms are shown underlined in 
the evolution equations (4) and (Al)), we use upwind-off-centered molecules (still 5 point 
4th order), with the direction of the upwinding chosen based on the sign of the radial shift 
vector (3 = (3^ . Tables II and III show this in detail. 

2. Grid Boundaries 

Our finite differencing near the grid boundaries may be described in either of two com- 
pletely equivalent ways, each useful in different contexts: 

Off-centered molecules: In this viewpoint, we use off-centered molecules for partial deriva- 
tives near the grid boundaries. To approximate a given w-coordinate partial deriva- 
tive (either 8^ or dww) near one of the boundaries, we choose from within one of the 
molecule families described below, the molecule with the closest possible offset (i.e. the 
closest possible centering or off-centering) to that used in the grid interior, subject to 
the constraint that the molecule not require data from outside the boundary. For a 
1st derivative d^ we use one of the off-centered 5 point 4th order molecules shown in 
part (a) of table III. For a 2nd derivative dwui we use either one of the off-centered 

5 point 3rd order molecules shown in part (b) of this table, or one of the off-centered 

6 point 4th order molecules shown in part (c). (We discuss the tradeoffs between these 
latter two choices in section VII B 3.) 

Extrapolation: In this viewpoint, we first logically extend the grid to include a few fictitious 
grid points beyond each boundary of the continuum problem domain, i.e. the fictitious 
points are at coordinates w = Wmin— 1, Wmin— 2, Wmin— 3, . . . inside the inner boundary, 
and w = Wmax+1, Wmax+2, Wmax+3, . . . outside the outer boundary. To approximate 
a given w-coordinate partial derivative {dw or d^w) near one of the boundaries, we 
then use the (assumed) smoothness of the grid function being finite differenced, to 
extrapolate values for the grid function at the fictitious grid points. For a 1st derivative 
dyj we use the 5 point 4th order (Lagrange) polynomial extrapolants^^ shown in part (a) 



^^We define the "order" of a Lagrange polynomial interpolant/extrapolant to be the order of the 
polynomial implicitly fitted to the known grid function values, the local truncation error of the 
extrapolated value being one order higher than this. Note that this definition differs from that 
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of table IV. For a 2nd derivative dww we use either those same extrapolants, or the 
6 point 5th order extrapolants shown in part (b) of this table. (Again, we discuss the 
tradeoffs between these latter two choices in section VII B 3.) Finally, using the newly- 
extrapolated values, we apply our grid-interior finite differencing scheme (including 
upwinding of Lie-derivative du, terms) at all the original (non-fictitious) grid points. 

By convolving the extrapolants in table IV with the various centered and off-centered 
molecules in table III, it's straightforward to show that the off-centered molecule and ex- 
trapolation viewpoints are in fact exactly equivalent - they yield identical finite difference 
equations, and hence identical results (modulo floating-point roundoff errors) . 

Our numerical code implements both viewpoints, and we find each to be useful in cer- 
tain circumstances: Off-centered molecules are most convenient when the actual molecule 
coefficients of a finite difference operator are desired. This case occurs, for example, when 
finite differencing the left hand side operators of elliptic PDEs such as those of coordinate 
conditions'*^ or the 3 + 1 initial data problem (cf. paper I). In contrast, extrapolation is most 
convenient when a finite difference operator need only be applied to a grid function, without 
necessarily explicitly computing any molecule coefficients. This case occurs, for example, 
when finite differencing the right hand side of the 3 + 1 evolution equations. 

3. 3rd versus 4th Order dyj^ Molecules Near the Boundaries 

In section VII B 2 we describe two slightly different finite differencing schemes, one where 
we use the same (4th) order of finite differencing for 2nd derivatives near the boundaries as 
elsewhere, and one where we use one order of accuracy lower (3rd order). Here we discuss 
the tradeoffs between these two schemes: 

Overall Accuracy: One would expect that when using only 3rd order accurate molecules for 
some terms in the 3 + 1 equations near the boundaries, the overall accuracy of our 
results would correspondingly be lowered to 3rd order. Near the outer boundary this 
is indeed the case. However, as discussed in section IX F, the dominant errors in our 
numerical solutions there are essentially continuum effects, not finite differencing ones, 
so lowering the finite differencing accuracy by one order makes little difference to the 
solutions' accuracies there. 



used for finite difference molecules, where the order is conventionally defined to be that of the local 
truncation error itself. 

'°The coordinate conditions described in section VI B can actually be solved independently at each 
grid point. However, in anticipation of future generalization to more general (elliptic) coordinate 
conditions, our numerical code uses a general-purpose elliptic solver to compute a and (3 at each 
time step. Each elliptic solution is done by first row scaling the resulting linear system to improve 
its numerical conditioning, then LU decomposing and solving it via LINPACK (Dongarra et al. 
(1979)) band matrix routines. 
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Near the inner boundary the situation is very different: By virtue of the black hole 
exclusion scheme, the light cones there point entirely into the black hole, so the future 
domain of dependence of a grid point near the inner boundary contains at most only a 
few other grid points. Gustafsson (1971, 1975, 1982) has shown theoretically, and Gary 
(1975a,b) has confirmed with numerical experiments, that for model problems with this 
"outflow boundary" ^^ causality, off-centered molecules 1 order of accuracy lower than 
the interior molecules may be used near the outflow boundary (the inner boundary 
in our system) without lowering the overall accuracy of the solutions. Unfortunately, 
these results hold only for simple model problems, not for nonlinear coupled hyperbolic- 
elliptic tensor PDEs such as the 3 + 1 evolution and coordinate equations. However, 
our numerical results in section IX E strongly suggest that these results do indeed 
generalize to the 3-1-1 equations: we find that the 3 + 1 field variables, and diagnostics 
computed from them, remain 4th order accurate even when using only 3rd order 
molecules for 2nd derivatives near the inner grid boundary. 

Efficiency and Implementation Convenience: In order to remain 4th order, the off-centered 
2nd derivative molecules in part (c) of table III must be one point larger (6 points) 
than the centered and 1st derivative molecules (5 points). When using the off-centered 
molecules viewpoint, we find that accomodating molecules of varying sizes generally 
adds only modest implementation complexities and inefficiencies. 

Correspondingly, the 5th order extrapolants in part (b) of table IV must be one point 
larger than the 4th order extrapolants in part (a). This extra size doesn't in itself 
cause any particular problem, but using different extrapolants for 1st and 2nd spatial 
derivatives^^ does have significant implications for the overall organization of a 3 + 1 
code: Because two different extrapolants may need to be applied to the same grid 
function(s), we can't just extrapolate all the to-be-finite-differenced grid functions once 
en masse at the start of each time step or evaluation of the 3 + 1 evolution equations' 
right hand side. Instead, we must perform a separate extrapolation for each individual 
finite differencing operation. ^^'^^ 



^^That is, "outflow" witfi respect to the computational domain. This terminology is common in 
the computational fluid dynamic literature, cf. footnote 6. 

^^Using the 6 point 5th order extrapolants of part (b) of table IV for 1st derivatives as well as 
2nd derivatives gives a violently unstable differencing scheme. 

^^This could be optimized to only re-extrapolate when changing from dw to dww or vice versa, 
but the basic issue remains unchanged. 

^^If we use the natural "pointwise" organization of the computations, where cache locality is 
optimized by interleaving algebraic and finite differencing operations so that we do as much com- 
putation as possible at one grid point before moving on to the next point, then the extrapolation 
is somewhat more awkward: For each individual finite differencing operation at each grid point, 
we must now test to see if we're close enough to the (a) boundary for extrapolation to be needed. 



24 



C. Time Integration 

Our time integration uses the method of hues. This technique is somewhat uncommon 
in numerical relativity, so we give a brief introduction to it in appendix B. 

As discussed in section VII A, for simplicity's sake we don't use any spatial adaptive 
gridding in this work. For the same reason, we use simple fixed-step fixed-order time in- 
tegrators. We have implemented both the classical 4th order Runge-Kutta method and 
a 4th order Adams-Bashforth-Moulton predictor-corrector method (Forsythe et al. (1977, 
section 6.4)). Each time step requires 4(2) right-hand-side-function evaluations for the 
Runge-Kutta (predictor-corrector) integrator respectively, i.e. 4 (2) complete evaluations of 
the coordinate conditions (14) and (15), the evolution equations (Al), and their respective 
boundary conditions (20), (21), and (17). 

We use a fixed time step, chosen to be slightly less than /T;(Ar)min, where k is an 0(1) 
parameter and (Ar)min = (Aw) / {dw / dr)^in is the r-coordinate grid spacing on the inner 
boundary (i.e. the smallest r-coordinate grid spacing anywhere in the grid).^^'^^ 

We have made test evolutions with a variety of k, values and grid resolutions to deter- 
mine our evolution scheme's empirical stability limit. ^^ For the spacetimes discussed in this 
paper, this limit is approximately k, = 1.5(0.5) for the Runge-Kutta (predictor-corrector) 
integrator. (This limit is weakly resolution-dependent in the sense of lower resolutions hav- 
ing a somewhat larger k limit, up to about a 20% increase at the lowest resolutions used 
here.) 

Comparing the two time integrators, the Runge-Kutta method is easier to program, 
more flexible (it allows easy changing of the time step), and when operating at the stability 
limit, about 50% more efficient. (That is, compared to the predictor-corrector integrator, 
the Runge-Kutta integrator requires twice as many right-hand-side-function evaluations per 



This adds a small extra overhead at each grid point, not just at the near-boundary points. 

^^In detail, we choose the time step At as follows: For compatability with the predictor-corrector 
integrator, we use the same (fixed) At for all time steps. (Changing At during an evolution is 
trivial for the Runge-Kutta integrator, but somewhat awkward for the predictor-corrector time 
integrator.) We don't do any time interpolation, so we (must) choose At to evenly divide all the 
coordinate times where field variables are to be output. To this end, for any real x > let [xj 
denote the largest number < x satisfying this latter condition. We then choose the time step to be 
At = Ik (Ar)mmJ . 

^^This is certainly not the only possible way to choose the time step. Another plausible scheme 
would be to normally use a time step k (Ar)min, except that if the next such step would take us 
past an output time, then take a single shorter Runge-Kutta step so as to just reach the output 
time. 

^'''Since the Einstein equations are nonlinear, the stability limit presumably also depends on the 
precise spacetime being evolved. In practice, though, this dependence doesn't seem to be very 
strong. 



25 



time step, but allows a maximum stable time step about 3 times larger.) All our results in 
this paper are from evolutions using the Runge-Kutta integrator with k, = 1.5. 

D. Comparison with Other Methods 

Our numerical methods differ in several ways from those usually used in 3 + 1 numerical 
relativity. Here we discuss some of the advantages and disadvantages of our key design 
choices: 

1. 4th Order Finite Differencing 

Given smooth grid functions and a reasonable grid resolution, 4th order finite differencing 
yields dramatically improved accuracy over the usual 2nd order, for (we find) only a modest 
increase in implementation effort and computational work per grid point. 

For time evolution problems, the most difficult part of 4th order finite differencing re- 
mains the same as that of 2nd order finite differencing: constructing (discovering) a stable 
differencing scheme. Contrary to some suggestions elsewhere (e.g. Press et al. (1992, sec- 
tion 19.3)), we do not find it any more difficult to construct stable 4th order differencing 
schemes than 2nd order ones. In fact, we find that allowing molecules to be larger than 
3 points (which is essential for any 4th order scheme except a nonlocal "compact" one 
(Ciment and Leventhal (1975); Hirsh (1975))) makes it easier to construct stable finite dif- 
ferencing schemes. This is because constructing a stable differencing scheme often requires 
considerable trial-and-error experimentation with various candidate schemes (e.g. Marsa and 
Choptuik (1996, section IV)), and larger molecules allow the formation of a wider variety of 
such candidates consistent with the necessary domains of dependence. (On the other hand, 
the requirement of 4th order accuracy disallows most averaging techniques and ADI-type 
schemes.) 

With 4th order finite differencing schemes of the type used here, there's no need to 
use a staggered grid of the type sometimes used with 2nd order finite differencing schemes 
(e.g. Choptuik (1986); Hawley and Evans (1989); Choptuik (1991)). We consider this a 
(modest) advantage, as staggered grids are somewhat awkward both in programming and 
in the eventual data analysis. 

It's also interesting to note that we haven't needed to use any type of causal differencing 
(Alcubierre (1993); Alcubierre and Schutz (1994)) to obtain a stable evolution scheme, even 
inside the horizon where the grid points are moving along spacelike world lines. We haven't 
investigated the reasons for this in detail, but we suspect it's due to the larger numerical 
domains of dependence inherent in our scheme's 5 and 6 point molecules. 

2. Richardson Extrapolation 

Another technique (which we haven't used here) for increasing the order of accuracy of 
any well-behaved (stable, convergent) finite differencing computation, is Richardson extrap- 
olation from lower order computations at multiple grid resolutions. (We briefly review this 
technique in appendix C.) This is a very powerful technique, and has been used successfully 
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by (e.g.) Choptuik et al. (1992) to obtain 4th order accuracy in spherically symmetric scalar 
field evolutions using 2nd order underlying finite differencing. 

For some simple model problems, it's in fact straightforward to show that such a Richard- 
son extrapolation scheme is in fact exactly equivalent (yields identical finite difference equa- 
tions) to a 4th order finite differencing scheme of the type we use here. For example, for the 
fiat-space linear scalar wave equation and 2nd order leapfrog finite differencing, Richardson 
extrapolating from results at a 2:1 ratio of grid spacings is precisely equivalent to the cen- 
tered 4th order molecules given in table III. However, it's not clear that this equivalence 
also holds for more complicated systems of PDEs like the 3 + 1 equations. 

Though still useful, Richardson extrapolation offers lesser benefits if parts of the differ- 
encing scheme are off-centered: As explained in appendix C, given a pair of computations 
at differing grid resolutions (say with a 2:1 ratio of resolutions), in general a single Richard- 
son extrapolation raises the order of accuracy by 2 if the finite differencing scheme is fully 
centered, but only by 1 if the scheme is even partially off-centered. Notably, in the presence 
of off-centering, Richardson-extrapolating a pair of 2nd order computations in general only 
yields a 3rd order result, or alternately three 2nd order computations with differing grid 
resolutions (say with a 4:2:1 ratio of resolutions) would be needed to obtain a 4th order 
result by Richardson extrapolation. 

Similarly, Richardson-extrapolating a 2:1 pair of our (partially off-centered) 4th order 
computations would yield "only" 5th order results. This might well be a useful scheme, but 
we haven't investigated it as yet. 

We should note one possible drawback to such a Richardson-extrapolation scheme: In 
practice, our numerical results in section IX D suggest that although our our 4th order 
evolution scheme yields qualitatively very accurate results at even modest grid resolutions, 
quite high grid resolutions are needed before quantitative 4th order convergence is fully 
attained, i.e. before the field variables accurately satisfy the convergence relationship (C2), 
or equivalently the Richardson error expansion (CI). 

For example, in figure 10, notice that the lOOolO — 200ol0 and 200ol0 — 400olO curves are 
still significantly different in shape for r ^ 70. This means that Richardson extrapolating the 
lOOolO and 200ol0 field variables would quite likely not yield significantly increased accuracy 
in this range of r; in fact it might well even decrease the accuracy there. The resolution 
would need to be increased by another factor of 2 (Richardson extrapolating the 200ol0 and 
400ol0 field variables) before we could be confident that Richardson extrapolation would 
work well (yield a significant accuracy improvement over the higher-resolution pair). 

3. Method of Lines 

We use the method of lines ( "MOL" ) to time-integrate the evolution equations. In com- 
parison to the more common "space and time together" (SATT) finite differencing schemes, 
MOL offers both advantages and disadvantages: 

Many SATT schemes in the numerical analysis literature are presented (only) for simple 
model problems, and are difficult to generalize to complicated and "messy" problems like the 
3 + 1 evolution equations. In particular, with SATT schemes it's often difficult to handle the 
2nd spatial derivatives on the right hand side of the 3 + 1 evolution equations, and schemes 
using substeps often encounter difficulty with obtaining suitable boundary conditions and/or 
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treating the coordinate equations. In contrast, MOL schemes usually generalize easily to 
"messy" problems. 

In the 3 + 1 context, MOL decouples the time integration from the detailed form of the 
evolution and coordinate equations. This decoupling is an advantage in that it isolates the 
complexity of the 3 + 1 equations in the spatial finite differencing, simplifying the program- 
ming, debugging, and stability analysis of a 3 + 1 code. On the other hand, this decoupling 
may also be somewhat of a disadvantage, in that MOL may give up an 0(1) factor in ef- 
ficiency relative to a equal-order SATT scheme. MOL also offers somewhat less flexibility 
to tailor the flnite differencing of individual terms in the 3 + 1 equations to their functional 
form for (say) improved stability and/or accuracy. 

MOL makes it relatively easy to use flnite differencing with higher than 2nd order accu- 
racy in both space and time. Although 4th order SATT schemes do exist (see, e.g., Kreiss 
and Oliger (1972, 1973); Oliger (1974); Ciment and Leventhal (1975); Hirsh (1975); Turkel 
et al. (1976); Kreiss (1978)), they generally don't apply to the 3 + 1 equations, due to the 
"messy problem" difficulties noted above. In contrast, constructing 4th order MOL schemes 
is straightforward, even for the full 3 + 1 equations. 

For almost all time evolution problems, and certainly for the 3 + 1 Einstein equations, 
the most difficult part of both SATT and MOL flnite differencing is the construction of 
stable differencing schemes. We flnd this to be about equally difficult for SATT and MOL; 
neither method seems to have a clear advantage here. 

We flnd (only) two signiflcant drawbacks to MOL: First, MOL schemes are less widely 
used than SATT schemes, and the MOL literature is still somewhat sparse. Second, the 
development of MOL adaptive gridding schemes analogous to the Berger-and-Oliger-type 
SATT adaptive gridding schemes noted in section VII A, remains an open research problem. 

VIII. DIAGNOSTICS 
A. The Energy Constraint 

One of our major diagnostics of the accuracy of our results is the numerically computed 
energy constraint C.^^ C has dimensions of (length)"^, but there are at least two different 
plausible ways to normalize it (make it dimensionless) : 

Some authors (e.g. Cook et al. (1998)) normalize C by taking the absolute value of each 
of its tensor- level terms, i.e. since C is given by (5a), they deflne 

Cabst = \R\ + \K,,K''\ + \K^\ + |167rp| (25) 

so that Creit = C/Cshst is a dimensionless measure of the extent to which the individual 
tensor-level terms in (5a) cancel out to yield (ideally) a zero sum C. 



^^C is a much more sensitive diagnostic for this purpose than the momentum constraint C*, 
because C depends on 2nd spatial derivatives of the 3-metric components, whereas C* only depends 
on 1st spatial derivatives. 
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However, each of these tensor-level terms is in fact computed as the sum of several 
individual terms involving partial derivatives of the state vector, and these "inner sums" 
may themselves involve nontrivial cancellations. For example, for an Eddington-Finkelstein 
slice of a mass-m Schwarzschild spacetime, at large r the computation of R by equation (B5) 
of paper I involves 2 orders in r of cancellations: R = 0{m'^ /r^) but some of the individual 
terms in that equation are 0{l/r'^). 

We can better measure all the cancellations involved in computing C by an alternative 
normalization, where we take the absolute value of each individual partial-derivative term. 
That is, by equations (B5) and (B15) of paper I, we have 
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so that Crcip = C jC^ii^ is a dimensionless measure of the extent to which the individual 
partial-derivative terms in (26) cancel out to yield (ideally) a zero sum C . 

As outlined above, Creit and Creip measure different things. For an Eddington-Finkelstein 
slice of a mass-m Schwarzschild spacetime, Cabst = 0{rr? jr'^^^ while Cabsp = ^(l/r^), so 
C'rcit/C'rcip = Oij-^ jrv?). That is, requiring Crdt to be below a given tolerance at large r is a 
much stricter criterion than requiring Crcip to be below that same tolerance. 

Since Crcip measures the extent of cancellation in the actual expression we use to com- 
pute C, we think it's a more realistic diagnostic of the numerical quality of our solutions. 
However, the definition of Creip depends on the details of our computational scheme: it would 
be difficult to unambiguously compare Crcip values for different schemes, or between (say) our 
spherically symmetric scheme and a generalization to axisymmetric or fully 3-dimensional 
spacetimes. In contrast, Crdt has a simple, elegant, and computational-scheme-independent 
definition. 

In describing our numerical results we generally give both diagnostics. Much of our 
discussion of them is common to both, so we often refer to Crci, meaning either Crdt or Crdp. 

B. Other Diagnostics 

Two other important diagnostic of the accuracy of our results are the changes in the 
numerically computed 3-metric and extrinsic curvature components when the grid resolu- 
tion is doubled in a convergence test. (We describe our convergence-test methodology in 
detail in appendix C) That is, defining Z[Ati;] to be the numerically computed value of 
the field variable Z using the grid spacing Aw, we use Agij = gij[Aw] — gij[Aw/2] and 
AKij = Kij [Aw] — Kij [Aw/2] as diagnostics of the accuracy of gij and Kij respectively. For 
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convenience we normally work with pointwise tensor "magnitude" norms of these diagnos- 
tics, {{Ag,,)) and {{AK,,)).^^ 

Our remaining diagnostics are mostly described in detail in paper I (section VIII and 
appendices B and C), so we only briefly summarize them here: To study the scalar field 
we use its 3-energy density p, its radial 3-energy and 3-momentum densities AnBp and 
4:71 Bj'^ respectively, the 4-Ricci scalar ^^^R = 87r(p — T) (which we usually normalize as 
^B ^'^^R = 4:7rB{p — T) to match AnBp and AirBj^), and the (Misner-Sharp) mass function m 
(Hayward (1994, 1996)).3°'3^ 

To study the black hole's growth and mass, we use the areal radius h of the apparent 
horizon(s) in each slice, and rribh, the mass function interpolated to the outermost apparent 
horizon's position. In a slight abuse of terminology we refer to rribh as the (time-dependent) 
"mass of the black hole" . As discussed in appendix B of paper I, in spherical symmetry any 
apparent horizon must have rarcai (= h in our coordinates) = 2m. For any diagnostic Z, Zh 
denotes the value of Z at (interpolated to) the or an apparent horizon position. 



IX. SAMPLE EVOLUTIONS 

We have made a number of test evolutions to investigate our code's stability, accuracy, 
and performance, and to study some of the spacetimes numerically generated by it. These 
evolutions include both vacuum (Schwarzschild) spacetimes, and dynamic spacetimes con- 
taining black holes surrounded by scalar field shells, but in the interests of brevity we only 
discuss the latter here. 

We discuss our initial data computation for these evolutions in detail in paper I. Briefly, 
we begin with an Eddington-Finkelstein slice of the unit-mass Schwarzschild spacetime, add 
a Gaussian to one of the field variable components (always P for the examples considered 
here), numerically solve the full 4- vector form of the York initial data algorithm to project 
the field variables back into the constraint hypersurface, and finally numerically coordinate- 
transform back to an areal radial coordinate. For the parameters used here, this algorithm 
results in a roughly Eddington-Finkelstein slice (with K generically nonzero and spatially 
variable everywhere in the slice) containing a black hole surrounded at some moderate radius 



^^Recall that for a symmetric rank 2 covariant tensor T^-, we define {{Tij)) = JTijT'^K 

^"^We actually use two different forms of the Misner-Sharp mass function: one computed directly 
from gij and Kij (we call this form mMSJ cf. appendix C of paper I), and the other computed 
by integrating the local mass density of Guven and O Murchadha (1995) (we call this form m^; 
cf. section VIII of paper I). As discussed in appendix B3 of paper I, at large r rriMS is quite sensitive 
to small numerical errors in gij and Kij, so we generally define m = m^, and refer to this as "the" 
mass function. 

^^When measuring the total mass of a slice, to reduce boundary-finite-differencing errors we 
actually measure the mass 3 grid points inside the outer grid boundary, i.e. we define mtotai = 

77l(w = Wmax-3). 
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by a roughly Gaussian scalar field shell. 

Table V summarizes various parameters of the evolutions discussed here. Notice that each 
model's name encodes the key numerical parameters (grid resolution and outer boundary 
position) and the basic physics (the input perturbation for the initial data solver). Most of 
our discussion of these models is is independent of their numerical parameters, so we often 
refer generically to the families of evolutions sharing common numerical parameters, and 
plot results for only one representative set of numerical parameters within each family. 

All the evolutions discussed here use an areal radial coordinate r, with the grid covering 
the range from w^iri = (rmin = 1-5) through one oi Wma.x = 4, 10, or 30 (r^ax ~ 248, 813, or 
2780). The numerical calculations use IEEE double (64 bit) precision for all floating-point 
arithmetic. 



A. Scalar Field Phenomenology 

We first consider the pqw5 evolutions, which show a black hole accreting a relatively 
thick and moderate-mass scalar field shell. On the initial slice the shell's radial density 
profile ATiBp has a standard deviation in r approximately 1.8 times the black hole radius, 
and the shell has approximately 0.66 times the black hole mass, i.e. 40% of the slice's total 
mass. 

Figure 3 shows the early-time evolution of the scalar field's radial 3-energy and 3- 
momentum densities AirBp and AttBj^, the normalized 4-Ricci scalar ^B ^^^R = 4:iTB{p — T), 
and the mass function m, for the 200.pqw5 evolution. (The other pqw5 evolutions don't 
differ noticeably at the scale and times of this figure.) 

It's clear that the initial slice's scalar field shell is actually a superposition of two 
momentarily-coincident shells, one propagating inward {j^ < 0), the other propagating out- 
ward [j^ > 0). The two shells have quite different masses, approximately 0.47 and 0.18 times 
the black hole mass (28% and 11% of the slice's total mass) respectively.^^ 

As the evolution proceeds, the outgoing scalar field shell propagates outward relatively 
intact, while the ingoing scalar field shell is partly captured by the black hole and partly 
scattered off the strong-field region's spacetime curvature. The scattered scalar field (visible 
in figure 3, as, for example, the positive pulse in AnBj^ at t = 15 near r = 20) then undergoes 
further scattering off the strong-field curvature. The net result of this "quasinormal-mode 
ringing" process is the formation of a sequence of successively-weaker outgoing scalar field 
shells following the original one. The outgoing scalar field shells all propagate out along 
asymptotically null geodesies, with asymptotically constant amplitudes in AiiBp and AttBj^, 
corresponding to ~ 1/r^ falloffs in p and j^. 

These features can be seen in figure 4, which shows the later-time evolution of the radial 



^^Given the nature of our initial data algorithm (cf. paper I) and the fact that the pqw5 models' 
initial perturbation is to P only, with Q remaining identically zero (cf. table V), the general form 
of a superposition of ingoing and outgoing shells is as expected. However, we don't know why the 
two shells' masses differ so much. 
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scalar field density and mass function for the 200ol0.pqw5 and 200.pqw5 models. ^^ (For 
the moment we consider only the 200ol0.pqw5 results; as discussed in section IX F, the 
200.pqw5 results are contaminated by outer-boundary errors.) For example (summarizing 
the scalar field shells' profiles by their approximate mean positions ± standard deviations 
in AnBp), at t = 200 the original outgoing scalar field shell is at r ^ 200 ± 10, while the 
quasinormal-mode-ringing shells are at r ^ 180 ± 20, 150 ± 20, 120 ± 20, and 80 ± 10. 

B. Quasinormal-Mode Ringing 

During the quasinormal-mode ringing process, at any fixed radius near the black hole, 
the field undergoes a sequence of damped (temporal) oscillations, whose frequency is theoret- 
ically predicted to be characteristic of the black hole mass. Figure 5 shows these oscillations 
for the 400ol0.pqw5 evolution. About 3 oscillation cycles are visible before the field settles 
down into a smooth decay at later times (discussed in detail in section IX G). Examining 
the ph minima/maxima times gives an oscillation period of about r = 32.5 ± 1, where the 
error estimate is dominated by the deviations of the oscillation from exact periodicity. 

From our results in section IX C, for this evolution the black hole mass is very nearly 
constant at m = 1.183 for t > 45 (where most of the oscillations occur), so our measured 
oscillation period corresponds to r = (27.5±l)mbh) in excellent agreement with the theo- 
retical prediction (Marsa (1995, section 4.4)) of r = 28.44 ?72bh- 

C. Black Hole Growth and Apparent Horizon Motion 

Because the black hole captures some of the scalar field, the black hole grows during the 
evolution. Figure 6 shows the time evolution of the horizon position h and the contained 
(black hole) mass m{h) for the 200.pqw5 model. Because the scalar field shell is relatively 
thick and of moderate mass compared to the black hole, the black hole grows relatively 
slowly and smoothly, and there's always only a single smoothly-moving apparent horizon. 

In contrast, consider the 400.pqwl model, in which the scalar field shell is relatively thin 
and massive. (On the initial slice, the shell's radial 3-energy density profile AnEp has a 
standard deviation in r of only about 0.43 times the black hole radius, and the shell has 
approximately 3 times the black hole mass, i.e. 75% of the slice's total mass.) 

The general phenomenology of this model's scalar field evolution is fairly similar to that 
of the 200.pqw5 model. However, as shown in figure 7, the time evolution of the apparent 
horizon position h and contained (black hole) mass m{h) is quite different. When the thin 
and massive scalar field shell falls through the horizon, the black hole grows very rapidly 
during a short time interval. In particular, during the time interval 19.13 ^ t ^ 19.71 there 
are three distinct apparent horizons present. Marsa (1995); Marsa and Choptuik (1996) 
have reported seeing similar behavior in their spherically-symmetric-scalar-field numerical 
evolutions, though with somewhat different initial data. As discussed in section IIIC, this 



^^Note that the sharp downward cusps visible in AiiBp in this figure do not represent zeros, just 
relatively narrow nonzero minima. 
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type of transient multiple-apparent-horizon behavior has important implications for how a 
black-hole-exclusion computational scheme should handle the inner boundary. 

In spherical symmetry, apparent horizons are (occur at positions given by) precisely the 
zeros of the horizon function H defined by (11). Generically, i/'s zeros (apparent horizons) 
always appear or disappear in pairs, corresponding to local minima or maxima in H just 
touching the r axis. This behavior is illustrated in figure 8, which shows H for the 400.pqwl 
evolution at times t = 19(0.1)20. Moreover, at the moment a local minima or maxima 
just touches the r axis, the H curve is generically tangent to the axis there, so at the 
moment of appearance or disappearance apparent horizons are generically moving with 
infinite coordinate speed dh/dt. (This is visible in figure 7 as the h{t) curve becoming 
vertical at the appearance and disappearance times t k, 19.13 and t ^ 19.71.) As discussed in 
section III C, this very fast movement of the apparent horizon also has important implications 
for how a black-hole-exclusion computational scheme should handle the inner boundary. 

It would be interesting to numerically compute the event horizon for this spacetime 
(presumably using the methods of Anninos et al. (1995a); Libson et al. (1996)) to compare 
it to the apparent horizons, but we haven't done this. 

D. Convergence Tests 

To assess the overall accuracy level of our numerical evolution scheme, we begin by 
considering the magnitude of the energy constraint, normalized as both Creit and Creip 
(cf. section VIII A). Figure 9 shows |Crcit| and |Crcip| in the inner part of the grid for 
the various-resolution ol0.pqw5 evolutions, at the time t = 100. (We discuss the outer part 
of the grid in section IX F.) 

Consider first the lOOolO.pqwS, 200ol0.pqw5, and 400ol0.pqw5 evolutions. For r ^ 130 
the Crei values are very small and are dominated by floating-point roundoff noise. ^^ At 
smaller radia (i.e. in part of the slice with nontrivial dynamics), the Crd values are still quite 
small: Even for our lowest resolution (the 100ol0.pqw5 evolution), Creit is still ^ 10^^ in mag- 
nitude throughout the grid, and the more realistically normalized Crdp (cf. section VIII A) 
is ^3 X 10~^. Moreover, we generally have excellent 4th order convergence of the |Crei| 
values: In particular, notice the very nice factor-of-16 ratio between the 200ol0.pqw5 and 
400ol0.pqw5 |Crci| values. 

However, the ICrei] values for the 800ol0.pqw5 evolution don't fully show the expected 
convergence: in regions where the scalar field is small but not negligible (e.g. r < 70), ICrdl 
(and the underlying \C\ from which (7^1 is computed) for this evolution shows significant non- 
smoothness (rapid variations from one grid point to the next), and is generally only a factor of 
~ 10 lower than for the 400ol0.pqw5 evolution. Although we don't have conclusive proof, we 
very strongly believe that these effects are both due entirely to floating-point roundoff errors. 
In support of this hypothesis, we note that the underlying \C\ values in question are very 
small (typically ^ 10~^^), and floating-point roundoff effects in the 3-1-1 equations generically 



^^RoundofF noise is easily recognizable by its extreme nonsmoothness (order-of-magnitude varia- 
tions from one grid point to the next). 
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grow as l/(Aw)^ as the grid is refined (since the equations contain spatial derivatives up 
to 2nd order). Furthermore, we note that the 800ol0.pqw5 evolution does show excellent 
4th order convergence near the peak of the outgoing scalar field shell (100 ^ r < 130), where 
finite differencing errors from the scalar field shell make \C\ somewhat larger (and hence we 
expect roundoff effects to be less severe in relative terms). 

In an attempt to better estimate the effects of roundoff errors on our numerical evolutions, 
we modified our numerical scheme to (optionally) simulate the effects of somewhat increasing 
these errors, by adding a tiny amount of white noise (the noise amplitude is 6 times the 
worst-case roundoff error of a single floating-point arithmetic operation) to the field variables 
several times per time step.^^ However, this turned out to not always offer reliable estimates 
of the roundoff effects'^^ (cf. our discussion in the next 2 paragraphs), and at least in its 
present incarnation, we don't recommend this technique to other researchers. 

The 800ol0n3.pqw5 evolution is identical to the 800ol0.pqw5 evolution except for hav- 
ing noise added in this manner, and figure 9 also shows this evolution's |Crei| values. As 
expected, the extra noise (which is uncorrelated from one grid point to the next) makes 
the 800ol0n3.pqw5 |Crci| values (and the underlying \C\ values) much less smooth than the 
800ol0.pqw5 values. This supports our hypothesis that the (lesser) non-smoothness of the 
latter \C\ values is due to "normal" floating-point roundoff effects. 

However, in some regions of the grid, the 800ol0n3.pqw5 |Crei| values (and the underlying 
\C\ values) are consistently smaller - by up to an order of magnitude - than the 800ol0.pqw5 
values! This is particularly prominent for 20 ^ r ^ 47 and rmin ^ r ^ 19, and is also 
noticable for 52 ^ r ^ 60. We don't understand why adding noise would make |Crci| (and 
|C|) smaller, though we suspect a dithering effect may be involved here. 

Another way to quantitatively assess the accuracy of our numerical evolution scheme is to 
consider 3-grid convergence tests (cf. appendix C) for the individual gij and Kij coordinate 
components. Figure 10 shows the results of some tests of this type, for the same olO.pqwS 
evolutions, again at the time t = 100. Here we have done a 3-grid convergence test on each 
gij and Kij component, but for conciseness we plot only the tensor magnitudes of the gtj 
and Kij differences between the different-resolution models, {(Agij)) and {(AKij)) .^''' 

The {{Agij)) and {{AKij)) values all show excellent 4th order convergence, except for 
some noise in the highest-resolution {{AKij)) values. This noise has the same characteristics 
described above for the 800ol0.pqw5 |Crci| values in figure 9, and from its qualitative form 
is almost certainly due to floating-point roundoff errors. As additional evidence for this, 
note that this noise increases just as we would expect when we add extra noise to the 
field variables, i.e. going from the ((i^jj[400ol0.pqw5] — i^jj[800ol0.pqw5])) values to the 



^^We describe the noise-addition algorithm in detail in appendix D. 

^^In hindsight this isn't unexpected: Actual roundoff errors are often highly discrete and/or highly 
correlated (Kahan (1996)), so (continuous and statistically independent) uniform white noise may 
not model them well, and in fact may systematically misestimate their actual effects (see also 
Kahan (1999)). 



^^Recall that for a symmetric rank 2 covariant tensor Tij, we define {{Tij)) = JTijT'^^ 
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((-ft'jj[400ol0.pqw5] — -ft'jj[800ol0n3.pqw5])) ones. There's no sign here of anything hke the 
anomalous behavior seen in the |Crei| resuhs. 



E. The Effects of 3rd versus 4th Order dww Molecules Near the Inner Boundary 

Our spatial finite differencing is generally 4th order (cf. section VII B), i.e. we approx- 
imate spatial partial derivatives by 4th order finite difference molecules. In particular, we 
approximate spatial 2nd partial derivatives by the 4th order molecules of part (c) of table III, 
or equivalently we use the quintic extrapolants of part (b) of table IV. 

However, as discussed in section VII B 3, there are some efficiency and implementation- 
convenience advantages to approximating spatial 2nd partial derivatives near the grid bound- 
aries by the 3rd order finite difference molecules of part (b) of table III, or equivalently the 
quartic extrapolants of part (a) of table IV. The 100e4.pqw5, 200e4.pqw5, and 400e4.pqw5 
evolutions use this alternative finite differencing scheme, but are otherwise identical to the 
lOO.pqwS, 200.pqw5, and 400.pqw5 evolutions (respectively). 

In general, the results of the e4.pqw5 evolutions are very similar to those of the cor- 
responding pqw5 evolutions. Compared to the pqw5 evolutions, the e4.pqw5 evolutions' 
errors (measured either by the constraints, or by 3-grid convergence tests on the gij and Kij 
components) are substantially larger at the innermost few grid points on the initial slice, but 
these larger errors damp out within a few time steps, and thereafter the e4.pqw5 evolutions 
have similar accuracies to the pqw5 evolutions. For example, figure 11 shows |Crcip| for the 
innermost part of the grid for the 200.pqw5 and 200e4.pqw5 evolutions. Notice that on the 
initial slice (t = 0), at the innermost grid point the 200e4.pqw5 evolution's error is about 
3 orders of magnitude larger than that of the 200.pqw5 evolution. However, even by t = 0.01 
the errors are within an order of magnitude of each other, and by t = 0.1 (3 normal-sized 
time steps for this evolution), or at larger radia, they're essentially identical. 

More generally, we find (details omitted for brevity) that the e4.pqw5 evolutions show full 
4th order convergence of |Crci| and other accuracy diagnostics. This extends the the model- 
problem analytical results of Gustafsson (1971, 1975, 1982) and the numerical experiments 
of Gary (1975a,b), to the full 3 + 1 initial data and evolution equations. 

F. Outer-Boundary Errors 

As discussed in section VIC, our computational domain only extends out to a finite 
outer boundary radius r = rmax- Conceptually, there are three distinct ways in which this 
introduces errors in our solutions: 

• Because our (continuum) outer boundary conditions only approximate the true r > 
^max physics of the 3 + 1 equations, "true" infinite-domain solutions of the equations 
won't exactly satisfy our (continuum) boundary conditions at r = rmax- In a slight 
abuse of language, we refer to this effect as the "inaccuracy" of our (continuum) outer 
boundary conditions. 

• Finite differencing errors occur when implementing our (continuum) outer boundary 
conditions. 
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• Floating point roundoff errors occur when implementing the finite differenced outer 
boundary conditions. 

For the spacetimes considered in this paper, the first of these processes - the inaccuracy of 
our outer boundary conditions - is always the dominant error source. 

Independent of what causes them, it's useful to (further) subdivide the outer boundary 
errors into "static mismatching errors" occurring at all times, and "dynamic reflection errors" 
occurring only when propagating waves in the field variables (e.g. those accompanying an 
outward-moving scalar field shell) reach the outer boundary. 

Static mismatching errors are caused by the falloff rates of the field variables at large 
r not precisely matching those we assumed in deriving our our outer boundary conditions 
(cf. section VIC and table I). The result of these mismatches is that our outer boundary 
conditions slightly perturb the field variables at the outermost grid point. In general this 
perturbation also violates the constraints. In a time evolution these errors spread throughout 
their domain of dependence, so the net result of static mismatching errors is a wave of 
(small) field variable errors and constraint violation, originating at the outer boundary and 
propagating inward at the speed of light. 

Figure 12 illustrates static mismatching errors. This shows the magnitude |Crci| of the 
energy constraint in the outer part of the grid for the 200ol0.pqw5 evolution, for times 
t = 0(50)200. The inward-propagating constraint- violation wave is clearly visible above 
the floating-point roundoff error "noise floor" . The constraint- violation wave has amplitude 
?«0.3 in |Crcit|, or ^3 X 10"^ in |Creip|. 

As a continuum effect, static mismatching errors are essentially independent of the finite 
difference grid resolution, ^^ and can only be made smaller by using more sophisticated outer 
boundary conditions or by moving the outer boundary radius r = rmax farther out. The 
effects of varying rmax can be seen by comparing the constraint-violation-wave amplitudes 
in (for example) the 200.pqw5, 200ol0.pqw5, and 200o30.pqw5 evolutions. Comparing these 
(details omitted for brevity), for this data the constraint violation near the outer boundary 
scales approximately as |Crcit| ~ (^max)'^'^^^'^'°^ or |Creip| ~ (rmax)"^'^^^"'^^- There's a slight 
Gibbs-type oscillation at the crest of the inward-propagating constraint-violation wavefront, 
but this doesn't grow significantly with time, and there's no sign of any finite differencing 
instability. The constraint violation stays approximately constant in amplitude in both 
|Creit| and |Creip| as it propagates inward. 

When outgoing waves in the field variables reach the outer boundary, then dynamic 
reflection errors also occur: one can think of these as being due to the outward-propagating 
waves not precisely matching our simple Sommerfeld-type outgoing-radiation condition (16). 
This mismatching causes some of the outgoing waves in both the geometry and matter fields 
to effectively reflect off the outer boundary and propagate back inwards. 

Dynamic reflection errors can be seen in figure 4 as the difference in AnBp between the 
200ol0.pqw5 and 200.pqw5 evolutions at times t > 350: The 200ol0.pqw5 evolution's outer 
boundary (at r^ax ~ 813) is far enough out so the matter hasn't yet reached it at these 



^^We have verified this (details omitted for brevity) by comparing the errors in the different- 
resolution ol0.pqw5 evolutions. 
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times, while the 200.pqw5 evolution's outer boundary (at rmax ~ 248) is relatively close 
in, so the matter encounters it relatively early in the evolution (at t ^ 225), resulting in 
the spurious inward-propagating scalar field shell visible at (for example)^^ r ^ 165 ± 15 at 
t = 350, r ^ 115 ± 15 at t = 400, and r ^ 65 ± 15 at t = 450. Comparing the amplitudes 
of the outgoing and refiected scalar field shells for the 200.pqw5 evolution, it's clear that 
only a very small fraction of the outgoing scalar field is refiected: the outer boundary's 
refiection coefficient (defined as the ratio of the refiected to the incident peak amplitudes), 
is approximately 10~^ in AiiBp, or 10~^ in the scalar field amplitudes P and Q. 

Like static mismatching errors, dynamic refiection errors are essentially a continuum 
effect, not a finite differencing artifact, and can only be made smaller by using more sophis- 
ticated outer boundary conditions or by moving the outer boundary radius r = rmax farther 
out. The effects of increasing rmax can again be seen by comparing the outer-boundary 
reflection coefficient in (for example) the 200.pqw5, 200ol0.pqw5, and 200o30.pqw5 evolu- 
tions. Comparing these (details again omitted for brevity), the outer-boundary reflection 
coefficient scales approximately as ~ l/(^max)'^ in 4:nBp, or ~ l/(rinax)^ in P or Q. 

The dynamic refiections in the geometry variables (as measured, say, by |Crci|) show 
quite different behavior: in these variables the reflection coefficient is approximately unity, 
i.e. the |Crci| values associated with the ingoing reflected matter are approximately the same 
as those associated with the outgoing wave prior to the reflection. Fortunately, these values 
are generally very small, so we don't consider this to be a serious problem. 

G. Late-Time Decay of the Scalar Field near the Black Hole 

After the quasinormal-mode ringing dies out, the scalar field near the black hole continues 
to decay. In the limit of very late times, both analytical perturbation theory (e.g. Nollert 
and Schmidt (1994); Gundlach et al. (1994a); Gomez et al. (1994)) and explicit numerical 
calculations (e.g. Gundlach et al. (1994b); Marsa (1995); Marsa and Choptuik (1996)) show 
that the scalar field at either any fixed coordinate radius r, or at the horizon position h{t), 
should decay as ~ 1/t^. In our computational scheme itself isn't readily accessible, but 
it's clear that P and Q should both display this same decay rate, and ^^^R and AnBp should 
decay at twice this rate, i.e. as ~ 1/t®. 

Figure 13(a) shows the late-time behavior of '^^^Rh and AnBph for the 400o30.pqw5 evo- 
lution. Both diagnostics decay roughly as ~ 1/t^^ for 200 ^ t ^ 500, then the decay rate 
steepens somewhat through t ^ 700. Both diagnostics then reach reach (nonzero) minima at 
t ~ 750, increase again to maxima at t ^ 950, and thereafter decay at later times, eventually 
asymptoting roughly to ~ 1/t^'^. (As discussed below, the sudden jump visible at t = 5000 is 
due to dynamic refiections off the outer boundary.) A closer examination shows that there 



^^Here we're again summarizing the shell's profiles by their approximate mean positions it 
standard deviations in AirBp), as in section IX A. 

Within the overlap region t < 500, the AirBph and ^'^'Rh values are identical between the two 
evolutions to within < 20 parts per million. 
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is also considerable further structure in these diagnostics: as shown in figure 13(b), they 
both vary systematically by factors of 2 or so in addition to the average decay rates just 
described. In other words, the diagnostics' actual decay rates aren't constant. 

This behavior clearly disagrees with both the perturbation-theory predictions, and pre- 
vious numerical results. We have not been able to identify a numerical artifact or error 
source which would account for any of the behavior we see: 

• Refiections from the outer boundary don't arrive back at the horizon until t ^ 5000,^^ 
so they shouldn't be relevant here. 

• Finite differencing errors don't appear to be responsible: A convergence test between 
the 100o30.pqw5, 200o30.pqw5, and 400o30.pqw5 evolutions shows that finite differ- 
encing errors in the 200o30.pqw5 AnBph and |*^^-'-R/i| values are generally ^0.1% for 
t < 2500,*^^ rising to (only) k, 0.3% by t = 5000, and the 400o30.pqw5 errors should 
be another factor of 16 smaller still. 

• Floating-point roundoff errors don't appear to be responsible: Comparing the 
200o30.pqw5 and 200o30n3.pqw5 evolutions shows that adding extra noise to the field 
variables changes both diagnostics by only ^ 1% for t < 3000, rising to (only) k, 10% by 
t = 4500. However, as noted in section IX D, our noise-addition technique may some- 



^^That is, since our initial data has compact support, with only negligible scalar field for r S> 20 
(e.g. AirBp < 10^^*^ for r > 75), we don't expect significant dynamic reflections from the outer 
boundary until the outgoing matter propagates out to rmax, and then the reflections must propagate 
back in to the horizon, for a total of t ^ 2rmax + 0(ln r) ^ 5500 ^^ until reflection effects should be 
visible at the horizon. It's thus somewhat peculiar that the reflections are already visible at t = 5000 
in figure 13(a). From examination of the detailed field variable profiles, the explanation for this 
is that as the reflected matter propagates back inwards, a series of low-amplitude precursor waves 
develops somewhat ahead of the main ingoing matter shell;^^ the t = 5000 point in figure 13(a) is 
caused by the arrival of the first of these precursor waves at the horizon. 

^^The O(lnr) term corrects for the fact that at large r, the outgoing speed of light is c+ = 
1 — 0(l/r) in our (asymptotically Eddington-Finkelstein) coordinates. 

^^The development of these precursor waves does seem to represent a type of instability in our 
evolution scheme, but we don't consider this to be a significant problem, because the waves are 
very weak (their peak amplitude in AirBp is only ~ 10~^^ at t = 5000), they seem to grow only for 
inward-propagating matter shells, and even then their growth is quite slow (their doubling time is 
roughly 500m). 

^^The convergence test does show much larger errors in the narrow time range 700 ^ t ^ 800, but 
even there the peak 200o30.pqw5 — 400o30.pqw5 difference is only ~ 1.5%. We haven't investigated 
the reasons for this narrow time range having an order of magnitude larger finite differencing error 
than other times, but it's interesting to note that this is the same time range as the minima in 
both the matter diagnostics themselves. 
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times systematically misestimate the effects of floating-point roundoff errors (Kahan 
(1996, 1999)), so we can't be certain these aren't significant error sources. 

Even though we haven't identified any cause for these results, we don't think they repre- 
sent genuine features of the true continuum physics. Rather, we suspect that they represent 
some as-yet-undetermined artifact in our numerical evolution scheme and/or our data anal- 
ysis techniques. 

X. CONCLUSIONS AND DIRECTIONS FOR FURTHER RESEARCH 

In this paper we present a number of results on the 3 + 1 numerical evolution of dynamic 
black hole spacetimes, using the "black hole exclusion" technique to avoid singularities. 

We first discuss the various design tradeoffs for the inner boundary placement in a black- 
hole-exclusion computational scheme: We argue that the (an) apparent horizon provides 
a natural reference point for placing the inner boundary, and we discuss the causality of 
such black-hole-excluded evolutions. We identify two main options for the inner boundary 
location: it may be coincident with the (an) apparent horizon, or it may be somewhat inside. 
We discuss the advantages and disadvantages of both options, concluding that at present 
it's not clear which is preferable. 

We then discuss the further complications which arise when our slices may contain multi- 
ple apparent horizons: In practice, these are often present only for a limited time interval in 
an evolution, and in particular the apparent horizon which a black-hole-exclusion evolution 
scheme has been using as a reference for its inner boundary placement, may at some point 
disappear. Assuming the slicing still contains another apparent horizons or horizons con- 
taining the disappearing one, we discuss how the evolution scheme might transition its inner 
boundary reference to the new "farther-out" apparent horizon. We again outline several 
possible design options, and conclude that further research is needed to investigate these 
and other possible ways of handling the disappearing-apparent-horizon problem. 

We then present a set of general criteria for choosing the coordinates in a black-hole- 
exclusion computational scheme. Like a number of other researchers (e.g. Marsa (1995); 
Marsa and Choptuik (1996)), we argue that Eddington-Finkelstein coordinates and their 
generalizations are ideally suited for use with the black-hole-exclusion technique. Unfortu- 
nately, although it's easy to generalize these coordinates for dynamic spherically symmetric 
spacetimes, they have no geometrically-preferred generalization to generic multi-black-hole 
spacetimes with no Killing vectors. 

We then turn to a detailed presentation of our black-hole-exclusion computational 
scheme itself. To summarize our scheme's key properties: 

• We require that a black hole be present in all our slices, including the initial data slice. 
As discussed in detail in paper I, we use the usual York algorithm to construct our 
initial data, solving the full 4-vector form of the York equations numerically because 
our slices have K nonzero and spatially variable. 

• Our formulation of the 3 + 1 equations uses a free evolution, where the constraints 
are not used explicitly in the evolution, but serve solely as diagnostics of the code's 
accuracy. We discuss two different ways to normalize the constraints. 
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• We place the inner boundary somewhat inside the apparent horizon (typically at ^ 75% 
of the initial slice's horizon radius), rather than on the horizon as is sometimes done. 
At present we maintain the inner boundary at a fixed spatial coordinate position 
throughout the evolution, though we plan to generalize this in the future. 

• We have experimented with a number of different slicing and spatial coordinate condi- 
tions, but in this paper we discuss only a single choice for these, the combination of a 
generalized Eddington-Finkelstein slicing and a constant-area radial coordinate. (Our 
initial data is such that this radial coordinate is in fact an areal one.) 

• We use finite differencing to discretize the 3 + 1 equations. However, the details of 
our finite differencing scheme are quite different from those commonly used in 3 + 1 
computations: 

— For spatial finite differencing, in the grid interior we use centered 5 point 4th order 
molecules, except for upwind-off-centered 5 point 4th order molecules for the 
Lie derivative (shift-vector advection) terms in the 3 + 1 evolution equations. 
Near the grid boundaries we use off-centered 5(6) point 4th order molecules for 
1st (2nd) derivatives. 

— We have also experimented with an alternate finite differencing scheme, in which 
we use 5 point 3rd order molecules for 2nd spatial derivatives near the grid bound- 
aries. This has some performance and implementation-convenience advantages, 
and still gives overall 4th order accuracy for the numerical solutions. 

— Contrary to some suggestions elsewhere (e.g. Press et al. (1992, section 19.3)), we 
do not find it any more difficult to construct stable 4th order differencing schemes 
than 2nd order ones. In fact, we find that in some ways it's easier, since allowing 
molecules to be larger than 3 points opens up a wider variety of possibilities for 
different candidate finite differencing schemes (cf. Marsa and Choptuik (1996, 
section IV)) consistent with the necessary domains of dependence. 

— Our time integration uses the method of lines, with a fully-explicit fixed-time-step 
4th order Runge-Kutta time integrator. 

— We use a single smoothly-graded nonuniform spatial grid, with the grid points at 
fixed (time-independent) coordinate positions. Our nonuniform gridding quali- 
tatively resembles using a logarithmic radial coordinates, but offers finer control 
over how the grid resolution varies with position; normally we adjust the param- 
eters so the grid has relatively high resolution in the strong-field region, dropping 
to a moderate (asymptotically constant) resolution at large r. 

— We don't use any type of staggered grid. 

— We don't use any type of "causal differencing" . 

Our numerical results show that this evolution scheme is stable and can run "forever": 
we have evolved for over t = 4000m with no difficulties. 

Even at the lowest resolution presented here {Aw = 1/100, yielding Ar/r ^ 3% at r = 
Sm and 2% at r = 20m), our evolution scheme is still very accurate. For example, after 
evolving the 100ol0.pqw5 initial data to t = 100 with nontrivial dynamics, ICrcitl ^ 10~^ 
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throughout the grid, and the more reahstically normahzed |Creip| (cf. section VIII A) is 
^3 X 10~^. The actual errors in the Qij and Kij components themselves, as determined by 
3-grid convergence tests, are also very small: the errors in Qij (Kij) are ^ 10~^ (3 x 10"'^) in 
most of the grid, rising to ~3 x 10^^ (3 x 10~^) near the horizon, and ~ 10"*^ (10 '^) ^^ the 
inner grid boundary. 

As expected from our scheme's 4th order accuracy, at higher resolutions all these error 
measures decrease very rapidly: Each doubling of the grid resolution reduces the errors by 
a factor of 16. We do see some floating-point-roundoff-error effects at very high resolutions 
{Aw = 1/800, yielding Ar/r f» 0.4% at r = Sm and 0.25% at r = 20m). For the resolutions 
we've used, roundoff effects are never serious, generally just slight increases and/or noisiness 
in \C\ and the 3-grid convergence tests' Kij error estimates. ^^ 

We use a simple pseudo-scalar Sommerfeld-type outgoing radiation outer boundary con- 
dition for the 3-1-1 evolution equations. This works quite well for the outgoing matter, 
probably because our nonuniform spatial grid allows us to uwe quite large outer bound- 
ary radia. For example, even at our innermost outer boundary position, rmax ~ 248, only 
about 10"'^ of the outgoing scalar field amplitude (i.e. about 10~® of the outgoing energy 
density) incident on the boundary is reflected back inwards, and this reflection coefficient 
falls approximately as ~ l/(^max)*^ (l/(^max)^) as the outer boundary is moved farther out. 

We have made a number of test evolutions to investigate our code's stability, accuracy, 
and performance. These evolutions include both vacuum (Schwarzschild) spacetimes, and 
dynamic spacetimes containing black holes surrounded by scalar field shells. Our results 
from these tests are generally consistent with theoretical predictions and past numerical 
studies, but there are a few surprises: 

When a relatively thin and massive scalar shell accretes into the black hole, for a short 
time 3 distinct apparent horizons are present. We present general arguments, and show 
a numerical example demonstrating, that apparent horizons should generically appear and 
disappear in pairs, and at their moment of appearance/disappearance they generically move 
with infinite coordinate speed (i.e. horizontally in the usual spacetime diagram). Such highly 
superluminal motion might cause difficulties for a black-hole-exclusion evolution scheme 
which attempts track the apparent horizon's location and use it as a reference point for the 
scheme's inner boundary placement. 

Although our measurements of quasinormal-mode ringing in our numerical evolutions 
show excellent agreement with both black hole perturbation theory and previous numerical 
calculations, our measurements of the late-time power-law tail of the field near the horizon 
strongly disagree with both perturbation theory and previous numerical results. We think 
this discrepancy is due to a numerical artifact of some sort in our computational scheme 
and/or data analysis techniques, but we don't have any specific explanation for it at this 
time. 

Our numerical code implementing the computational scheme described here and in pa- 



^^In one test evolution we do see a very peculiar effect where artificially adding additional white 
noise to the field variables several times per time step throughout an evolution, causes \C\ to 
decrease by up to an order of magnitude in some parts of the grid! We suspect this may be a 
dithering effect, but we don't understand the details of this. 
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per I is available from the author's web site.^^ The code is written in ISO C (about 35K lines, 
although this includes considerable code for additional computations not described either 
here or in paper I), and should be easily portable to any modern computer system. It may 
be freely modified and/or redistributed under the terms of the GNU General Public License. 
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APPENDIX A: EVOLUTION EQUATIONS FOR THE SPHERICALLY 
SYMMETRIC SCALAR FIELD SYSTEM 

We have tabulated the detailed 3 + 1 equations for the constraints and all our diagnostics, 
written explicitly in terms of our state variables A, B, X, Y, P, and Q, in appendix B of 
paper I and section VIII A of this paper. Here we tabulate the corresponding equations for 
the time evolution of the state variables: 



dtA= -2aX + (3drA + 2{drp)A (Ala 

dtB= - 2aY + /3drB 

dfX = — drr-a + \{dra) 



dtB= - 2aY + l3drB (Alb) 

drA 




, idrA){drB) XY X^ 

. „_. , , , nOi : \- 2a a— r 

B ^ \ B J ^ AB B A 

+ (3dr^ + 2{dr0X 

- 2aP^ (Ale) 



d,Y= -Udr-a)^ 



-I drrB -I {drA){drB) XY 

- 2"^4- + 4" ^^ + " + "^4- 

+ f3 drY (Aid) 
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dtP = {dra)Q + adrQ 

+ {drP)P + p drP (Ale) 

AT 1^ 

+ a—Q + 2a—Q 

A JD 

+ PdrQ, (Alf) 

where the partial derivatives in the shift vector Lie derivative terms are underhned for 
reasons discussed in section VII B. 



APPENDIX B: A BRIEF INTRODUCTION TO THE METHOD OF LINES 

In this appendix we give a brief introduction to the method of lines ( "MOL" ) for parabolic 
or hyperbolic PDEs. This technique has a number of desirable features, but it's not widely 
known and the literature on it is quite sparse. 

The best general MOL references we're aware of are Hyman (1976, 1979). Other very 
good MOL references are Madsen and Sincovec (1973, 1974); Vemuri and Karplus (1981); 
Hyman (1989). Unfortunately, all of these except Vemuri and Karplus (1981) are some- 
what inaccessible, and none go into the level of detail necessary for a full understanding 
of the design choices in developing a method of lines code. We have discussed MOL at 
this latter level of detail in Thornburg (1993, sections 7.3.2-7.3.9). Other useful MOL ref- 
erences include Liskovets (1965) (reviews much of the earlier Soviet work on this topic). 
Carver (1976), Schiesser (1991) (very elementary), and Gustafsson (1971, 1975, 1982); Gary 
(1975a,b) (technical discussions of finite differencing techniques for boundary conditions). 

The basic idea of MOL is simple: We initially finite difference^^ only the spatial deriva- 
tives in the PDE, keeping the time derivatives continuous. This yields a set of coupled 
ODEs'*^ for the time dependence of the field variables at the spatial grid points. (In the 
terminology of relativity, these ODEs give the time dependence of the field variables along 
the spatial-grid-point world lines.) A suitable ODE integrator is then used to time-integrate 
these ODEs. 

A simple example should help to clarify this: Consider the flat-space linear scalar advec- 
tion equation written in 1st order form, 

dtu = —cdxU (Bl) 

on the domain (t,x) G [0, oo) x [0, 1), with the (smooth) coefficient function c being every- 
where positive (so that propagation is solely rightward), subject to the periodic boundary 
condition 



'^^The method of lines can also be used with finite element or pseudospectral discretizations, but 
for simplicity we only discuss finite difference methods here. 

The coupHng is due to the spatial finite differencing; the example below (starting with (Bl)) 
should clarify this. 
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u{t, x=0) = u{t, x=l) . (B2) 

To treat this problem by MOL, we first discretize the spatial dimension with the uniform 
grid Xk = k Ax for A; = 0, 1, 2, ... , K—1, where Ax = 1/K for some integer K. (For the 
moment we keep the time dimension continuous.) Introducing the usual notation Uk = u{xk) 
and Cfc = c{xk), and approximating the spatial derivative in the evolution equation (Bl) with 
the usual centered 2nd order finite differencing molecule, we obtain the coupled system of 
ODEs 

r% '?^fc+l mod iiT "iffe-l mod iiT fc i n i o r^ i \ f-nn\ 

dtUk = -Ck -^ (for A; = 0, 1, 2, ... , K-1) (B3) 

for the time dependence of the state variables {uk}. These ODEs can then be integrated 
using any suitable ODE integrator. 

Although the time integration of ODE systems is generally most easily, accurately, and 
efficiently done using modern adaptive methods, ^^ for MOL ODE systems simpler methods 
are often appropriate. In particular, in MOL there's normally little point in making the 
time integration much more accurate than the spatial finite differencing. Assuming a PDE 
whose solutions vary roughly equally rapidly in space and time, then for the (common) case 
where no spatial adaptive gridding is done and the spatial finite differencing is of fixed order, 
this suggests using a correspondingly simple fixed-step fixed-order time integrator. Hyman 
(1989) and Thornburg (1993, section 7.3.9) discuss the choice of MOL time integrators in 
more detail. 

Despite its superficial differences, MOL is actually closely connected to the more common 
"space and time together" (SATT) finite difference methods for PDEs: Essentially all ODE 
integrators use finite differencing in the time dimension, so the overall result of applying 
such an integrator to a system of MOL ODEs, is to (implicitly) construct and solve a 
(complicated) system of finite difference equations in space and time. 

For example, if we were to use the leapfrog time integrator 

u(t + At) = u{t- At) + 2At{dtu){t)+0{{Atf) (B4) 

(where u denotes the state vector {uk}) to time integrate the MOL ODEs (B3), it's easy to 
see that at each time level t„ = n At we would obtain the finite difference equations 

^^ = -Ck ^^^ {iork = 0,l,2,...,K-l) (B5) 



■^^These are discussed in detail by (e.g.) Gear (1971); Shampine and Gordon (1975); Gear (1981); 
Gupta et al. (1985); Byrne and Hindmarsh (1987), and are implemented in numerical codes such 
as (e.g.) DE/STEP (Shampine and Gordon (1975)) ODEPACK/LSODE (Hindmarsh (1993)), 
and RKSUITE (Brankin et al. (1992)). However, such methods and codes don't yet handle the 
varying-size ODE systems which would result from applying MOL to a Berger-and-Oliger-type 
adaptive gridding scheme (cf. section VII A). More generally, the integration of such adaptive 
gridding schemes into MOL remains an open problem. 
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where we use the notation nUk = u(t=tn,x=Xk). These are precisely the equations for 
the standard 2nd order centered-in-time centered-in-space ("CTCS") leapfrog SATT finite 
differencing scheme, applied to the original PDE (Bl). 

Corresponding to the equivalence of a MOL scheme to a suitable SATT scheme, the 
usual CFL stability limit of an explicit SATT scheme has an exact analog in MOL: When 
time integrating a system of MOL ODEs with an explicit ODE integrator, the maximum 
time step will be limited by stability considerations (Thornburg (1993, section 7.3.9)) to 
precisely the same value as the equivalent SATT scheme's CFL limit. For example, the 
scalar-wave-equation finite difference equations (B5) have the stability limit |cAt| < Ax, 
which may be viewed either as the usual CFL limit for the CTCS leapfrog SATT scheme, 
or the ODE-integration stability limit for the leapfrog time integrator (B4) applied to the 
system of ODEs (B3). 

Now consider a variant form of our example problem, where we extend the spatial domain 
from [0, 1) to [0, 1], and introduce the left boundary condition 

u{t,x=0) = smt. (B6) 

(There is no right boundary condition: The system (Bl) and (B6) is well-posed without 
one, since no part of the problem domain causally depends on the right boundary.) 

To treat this variant by the MOL, we first extend the grid to include xk = 1- We 
then time-differentiate the left boundary condition (B6) and combine it with the interior 
evolution equation (Bl) to obtain the "merged" evolution equation 

r. I cost (if X = 0) /„_^ 

*"=|-c4„ (,fx>0)- (^^' 

Finally, we spatially finite difference this in the same manner as before for the grid 
interior, but using the off-centered finite difference molecule 



d^ 



1 



+ 1 -4 +3 



+ OiiAxf) (B8) 



2 Ax 
at the right boundary grid point, to obtain the MOL ODEs 

dfUo = cost (B9a) 

dtUk = -0^ """+^"^'-^ (for A; = 1, 2, 3, ... , K-l) (B9b) 

OtUK = -Ck TT^ • (B9c) 

I Ax 

These can again be integrated using any suitable ODE integrator. (Note that the leapfrog 
ODE integrator (B4) would not be suitable here, in fact it would be unconditionally unstable 
for this problem. However, this is "just" a problem with this time integrator; these ODEs 
can easily be stably integrated by using a time integration scheme with better stability 
properties. For example, any of the integrators discussed in section VII C would suffice 
here.) 
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APPENDIX C: CONVERGENCE TESTS AND RICHARDSON 

EXTRAPOLATION 

In this appendix we briefly review our methodology for convergence testing and also 
briefly discuss the related technique of Richardson extrapolation. We have previously dis- 
cussed some of this material in more detail in appendix E of paper I. 

1. Convergence Tests 

As has been forcefully emphasized by Choptuik (1986, 1991); Choptuik et al. (1992), a 
careful comparison of the errors in approximating the same physical system with the same 
algorithm, but at different grid resolutions, can yield a great deal of information about, and 
very stringent tests of, a computational scheme's numerical performance and correctness. 

In the present context, consider some diagnostic grid function Z, and consider first the 
case where its true (continuum) value Z* is known a priori. (For example, in section IX D 
we take Z to be \C\ at time t = 100 for the olO.pqwS evolutions, so Z* = 0.) Suppose we 
have a pair of numerical computations of Z, identical except for having a 2:1 ratio of grid 
spacings Aw. As discussed in detail by Choptuik (1991), if all the field variables are smooth 
and the code's numerical errors are dominated by truncation errors from nth order finite 
differencing, then the numerically computed Z values must satisfy the Richardson expansion 



Z[Aw] =Z* + {AwYf + 0((A«;)"+2) (Cla) 

Z[Aw/2] = Z* + {Aw/2Yf + 0((Aw;)"+2) (Clb) 

at each grid point, where Z[Aw] denotes the numerically computed Z using the grid spacing 
Aw, and / is an 0(1) smooth function depending on various high order derivatives of Z* 
and the field variables, but not on the grid resolution. [We're assuming centered finite 
differencing here in writing the higher order terms as 0((Aw)""'"^), otherwise they would 
only be 0{{Aw)'^^^)] Neglecting the higher order terms, i.e. in the limit of small Aw, we 
can eliminate / to obtain the "2-grid" convergence relationship 

Z[Aw/2]-Z* = ^{z[Aw]-Z*), (C2) 

which must hold at each grid point common to the two grids. 

If Z* isn't known ahead of time (for example, in section IX D we also take Z to be 
an individual coordinate component of gij or Kij at time t = 100 for the olO.pqwS evo- 
lutions), then we consider a triplet of numerical computations of Z, identical except for 
having a 4:2:1 ratio of grid spacings Aw. Analogously to the previous case, we now have 
the Richardson expansion 

Z[Aw] =Z* + {AwYf + 0((A«;)"+2) (C3a) 

Z[Aw/2\ = Z* + (Aw/2)"/ + 0((Aw)"+2) (C3b) 

Z[Aw/A] = Z* + (Aw/4)"/ + 0((Aw)"+2) (C3c) 
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at each grid point. Again neglecting the higher order terms, i.e. in the hmit of small Aw, 
we can eliminate both / and Z* to obtain the "3-grid" convergence relationship 

Z[Aw/2] - Z[Aw/A] = ^X^[Aw] - Z[Aw/2]) , (C4) 

which must hold at each grid point common to all three grids. 

To assess how well numerical data satisfy one of the convergence relationships (C2) 
or (C4), we plot the differences Z[Aw]-Z* and Z[Aw]-Z* (2-grid) or Z[Aw/2]- Z[Aw/A] 
and Z[Aw] — Z[Aw/2] (3-grid) as a function of position, on a common logarithmic scale. If 
(and only if) the data satisfy the convergence criterion, the two plots will be identical except 
for a vertical offset corresponding to a factor of 2". 

It's important to note that we apply both the 2-grid and 3-grid convergence tests point- 
wise, i.e. independently at each grid point common to the different-resolution grids. This 
makes these analyses much more sensitive than a simple comparison of gridwise norms, 
which would tend to "wash out" any convergence problems occurring at only a small subset 
of the grid points (e.g. near a boundary). 

Notice that for both the 2-grid and the 3-grid convergence tests, the parameter n, the 
order of the convergence, is (or at least should be) known in advance from the form of 
the finite differencing scheme. (For our computational scheme n = 4.) Thus the factor 
of 2" convergence ratio in (C2) or (C4) isn't fitted to the data points, but is rather an 
a priori prediction with no adjustable parameters. A convergence test of either type is thus 
a very strong test of the validity of the overall finite differencing scheme and the Richardson 
expansion (CI) or (C3). 

2. Richardson Extrapolation 

So far we have described how to use the Richardson error expansions (CI) and (C3) to 
assess the accuracy of our numerical results. Alternatively, we can use these same expan- 
sions to (in many circumstances) improve this accuracy, by the technique of "Richardson 
extrapolation" : Suppose we once again have a pair of numerical computations of Z, identical 
except for having a 2:1 ratio of grid spacings Aw, but now suppose further that we don't 
know the true (continuum) value Z*. If we assume that the Richardson error expansion (CI) 
does indeed hold, then we can solve directly for Z*: 

nn 

Z* = Z[Aw] - -(z[Aw] - Z[Aw/2]) + 0((A«;)"+') (C5a) 

= Z[Aw/2] - ^-L-^[z[Aw] - Z[Aw/2]) + 0((A«;)"+^) (C5b) 

This is a remarkable result: the computed Z* is 2 orders more accurate than any of the 
numerical computations used to calculate it! [If the differencing scheme isn't fully centered, 
then the computed Z* is only accurate to 0((A-U7)"+^).] 

Notice that when doing Richardson extrapolation, we must assume that the Richardson 
error expansion (CI) does indeed hold. However, in practice, although we can't verify this 
directly if Z* is unknown (the only case where Richardson extrapolation is interesting), we 
can verify it indirectly by doing a 2-grid convergence test on (say) the energy constraint C 
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for the same 2:1 pair of computations: if this shows good convergence at the expected order, 
then it's a fairly safe assumption that other diagnostics hke Z (from the same numerical 
computation) also show similar convergence. 

However, there is one drawback to Richardson extrapolation: there's no longer any way 
to estimate the remaining 0((Atf;)"'+^) error [0((Aw)"+^) for the non-centered case] in the 
Richardson-extrapolated Z* result. We discuss the use of Richardson extrapolation further 
in section VII D 2. 

APPENDIX D: ADDING NOISE TO THE EVOLUTION 

Although the effects of floating-point roundoff errors are fairly easy to model for each in- 
dividual floating-point arithmetic operation (Knuth (1997, section 4.2.2), Goldberg (1991)), 
it's not easy to estimate the net (cumulative) effects of these errors on an entire computa- 
tional scheme such as ours. 

In the spirit of convergence testing (cf. appendix C), the obvious way to try to estimate 
the effects of floating-point roundoff errors would be to repeat our computations with a higher 
floating-point precision, and see how much various diagnostics change. Unfortunately, this 
would be somewhat awkward to do, since we already use IEEE double (64 bit) precision for 
all floating-point arithmetic, and this is the highest precision generally available. Instead, 
to estimate roundoff effects we repeat our computations with a simulated lower precision. 
The obvious choice of IEEE single (32 bit) precision is much too low (results with this 
precision would be quite inaccurate), so we use a software simulation of a slight lowering of 
the precision, more precisely a slight increase in floating-point roundoff errors. 

In detail, we simulate increased floating-point rounding errors as follows: Each time our 
time stepping algorithm calls for a right-hand-side function evaluation (i.e. 4 times per time 
step, cf. section VII C), we add a tiny amount of uniformly distributed white noise^° to each 
component of the time-integration state vector (A, 5, X, F, P, and Q). We then use the 
added-noise state vector as usual in computing the Einstein equations and integrating them 
ahead to the next time level. Note that this means that the added noise is cumulative: after 
A^ time steps we have added a total of 4iV white-noise samples (interspersed throughout the 
time evolution) to each state-vector component at each grid point. 

Each time we add noise, we scale its amplitude independently at each grid point to that 
of the corresponding state-vector component there: We scale the noise to be uniform over 
the range ±3 ulp for that component, where a "ulp" denotes a "unit in the last place" 
(Knuth (1997, section 4.2.2), Goldberg (1991)), the value of the least signiflcant fraction 
bit of a binary floating-point number. For comparison, floating-point roundoff errors are 
often modelled as adding white noise uniform over the range ±i ulp for each floating- 



^°We generate the white noise via the Unix random (3) random number generator. This is a 
high-quality generator using a nonlinear additive feedback algorithm (Knuth (1997, section 3.2.2)); 
it does not suffer from the usual approximate-linear-dependence or low-order-bit-nonrandomness 
problems which afflict most linear-congruential random number generators (Knuth (1997, sec- 
tion 3.3.4)). 
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point arithmetic operation (Knuth (1997, section 4.2.2), Goldberg (1991)). [However, it's 
important to realise that such a model has serious flaws: actual floating-point roundoff 
errors are often highly discrete and/or highly correlated (Kahan (1996)), so naive statistical 
analyses based on modelling them as uncorrelated continuous random variables can be highly 
misleading (see also Kahan (1999)).] 
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FIG. 1. This figure shows the hght cones and the Eddington-Finkelstein slicing of Schwarzschild 
spacetime. Part (a) is plotted in Kruskal-Szekeres coordinates {u,v), and shows the r = 1.5(0.5)3.5 
surfaces^^ and the Eddington-Finkelstein t = —8(1)4 slices, with the t = —4(1)4 slices labelled. 
Part (b) is plotted in Eddington-Finkelstein coordinates (r, t), and shows the Eddington-Finkelstein 
t = —1(0.5)3 slices and the r = 1.5 and r = 2 (horizon) surfaces. The legend is common to both 
parts of the figure. 

FIG. 2. This figure shows the behavior of our nonuniform gridding coordinate w. Part (a) 
shows the actual w coordinate, while part (b) shows the grid spacing Ar for grids with resolutions 
of (from top to bottom) Aw = 0.01, 0.005, 0.0025, and 0.00125 (we refer to these as "100", "200", 
"400", and "800" evolutions respectively). The diagonal dashed lines labeled along the top and 
right of part (b) show lines of constant relative grid spacing Ar/r. In both parts of the figure the 
vertical dashed lines show the outer grid boundaries for li^max = 4, 10, and 30. 

FIG. 3. This figure shows the scalar field's radial 3-energy and 3-momentum densities inBp 
and AirBj"^ , the normalized 4-Ricci scalar -^B ^'^'R = AirB^p — T), and the mass function m, for the 
200.pqw5 evolution at times t = 0(5)50. The vertical dashed line near the left side of each subplot 
shows the apparent horizon position h. Notice that the left vertical scale changes by a factor of 3 
between the t = 0-25 and t = 25-50 subfigures; for comparison the t = 25 subfigure is shown 
twice, once at each scale. 

FIG. 4. This figure shows the scalar field's radial density profile AnBp (on a logarithmic scale) 
and the m function (on a linear scale) for the 200ol0.pqw5 and 200.pqw5 evolutions at times 
t = 0(50)500. Notice the change in scales for r and AirBp (but not m) compared with figure 3: 
the box in the t = and t = 50 subplots shows the area plotted in that figure. The 200.pqw5 
AnBp curves in the t < 300 subfigures, and the 200.pqw5 m curves for all the subfigures, are at 
this scale indistinguishable from, and hence hidden under, the 200ol0.pqw5 curves. The 200.pqw5 
iirBp curves in the t > 350 subplots show significant outer boundary refiections; we discuss these 
in section IX F. 

FIG. 5. This figure shows the quasinormal-mode oscillations of the 4-Ricci scalar at the horizon, 
^ 'Rh (upper curve), and the energy density at the horizon, ph (lower curve), for the 400ol0.pqw5 
evolution at times < t < 200. (For visual clarity the plot actually shows p^/lOO instead of ph-) 
Note that while ^^'Rh changes sign (shown by the different plot symbols) in the oscillations, ph is by 
definition always positive semi-definite, and is in fact strictly positive at all times in this evolution. 

FIG. 6. This figure shows the apparent horizon's position h and contained mass m{h) for the 
200.pqw5 evolution at times < t < 50. (At later times the black hole doesn't grow significantly.) 

FIG. 7. This figure shows the apparent horizon's position h and contained mass m(h) for the 
400.pqwl evolution, (a) for the time range < t < 50, and (b) at an expanded scale for the time 
range 19 < t < 20. (The uneven density of points is to ensure good resolution of the near-vertical 
parts of the curve in part (b).) 

FIG. 8. This figure shows the left hand side of the apparent horizon equation, the horizon 
function H defined by (11), for the 400.pqwl evolution at times t = 19(0.1)20. As discussed in the 
text (section IXC), apparent horizons correspond to zeros of H. 
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FIG. 9. This figure shows the energy constraint for the olO.pqwS evolutions, at time t = 100. 
Part (a) shows the energy constraint normahzed as Creit, while part (b) shows it normalized as 
Creip; the vertical scale is the same in both parts. For reference, both subfigures also show the 
scalar field's scaled radial density AirBp (which at this scale doesn't differ appreciably between 
the different-resolution evolutions), at the same time. The vertical dashed line near the left side 
of each subfigure shows the apparent horizon position h. The scale bars show a factor of 16 in 
|Crei|, for comparison with the vertical spacings between the energy constraint plots at the different 
resolutions. 



FIG. 10. This figure shows the tensor magnitudes of the differences in the 3-metric and ex- 
trinsic curvature, {{Agij)) and {{AKij)) (parts (a) and (b) respectively), between pairs of ol0.pqw5 
evolutions which are identical except for a factor of 2 difference in grid resolution, all at time 
t = 100. For reference, it also shows the scalar field's scaled radial density AirBp (which at this 
scale doesn't differ appreciably between the different-resolution evolutions), at the same time. The 
vertical dashed line near the left side of each subfigure shows the apparent horizon position h. The 
scale bars show a factor of 16 in {{Agij)) and {{AKij)), for comparison with the vertical spacings 
between the plots at the different resolution pairs. 



FIG. 11. This figure shows the normalized energy constraint |Creip| for the innermost part of 
the grid in the 200.pqw5 and 200e4.pqw5 evolutions, at times t = 0, 0.01, 0.1, 1, 10, and 100. The 
vertical dashed line near the right side of each subplot except the t = 100 one shows the apparent 
horizon position h; for t = 100 the horizon has expanded to include the entire subplot. 



FIG. 12. This figure shows the effects of static mismatching errors in the outer boundary 
conditions. It shows the energy constraint in the outer part of the grid at times t = 0(50)200, 
for the 200ol0.pqw5 evolution. Part (a) shows the energy constraint normalized as Crdt, while 
part (b) shows it normalized as Creip; the legend and the axis scales are common to both parts. 



FIG. 13. This figure shows the decay of the scalar field at the apparent horizon in the 
400o30.pqw5 and 400ol0.pqw5 evolutions, for times 150 < t < 5000. (At earlier times the field is 
dominated by quasinormal-mode ringing, cf. section IX B and figure 5.) Part (a) shows the energy 
density at the horizon, p/^, and the magnitude of the 4-Ricci scalar at the horizon, |(^)i?;i|, along 
with ~ 1/t^^ and ~ 1/t falloff lines. Part (b) shows the ratios of the same data points to the falloff 
lines plotted in part (a), together with similar ly-ratioed falloffs r-^t '"^ and ~t faster/slower 
to give an indication of the uncertainties in the nominal falloff rates. Part (a) shows data from 
the 400o30.pqw5 evolution; part (b) shows both this same data, and also the (much more densely 
sampled) t < 500 data from the 400ol0.pqw5 evolution. ^'^ 
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TABLE I. This table gives the parameters m and n used for the evolution outer boundary 
condition (17). 
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TABLE IL This table gives the offsets of the finite difference molecules we use in various parts 
of the grid. For the shift vector Lie derivative terms in the evolution equations (these terms are 
shown underlined in the evolution equations (4) and (Al)), the molecule offsets also depend on the 
sign of the shift vector. The molecules themselves are given in table IIL 
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TABLE III. This table gives the coefficients of all of our finite difference molecules as a function 
of their offsets, (a) for 5^,, (b) and (c) for dww (As discussed in the text, parts (b) and (c) 
correspond to using the extrapolants of parts (a) and (b) of table IV, respectively.) The i+k 
column headings denote grid points for molecules applied at grid point i. (The application point 
is also denoted by underlining.) In part (b), the centered molecule has an error of 0((Ati;)^), 
while the off-centered molecules have errors of ©((A-u;)^). In part (c), notice that the off-centered 
2nd derivative molecules must be 1 point larger than the centered one (6 points instad of 5) in 
order to still have 4th order local truncation error. 
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Part (a): 5 point 4th order (quartic) extrapolants: 

/^T1 = + 5/.- 10/,±i+ 10/,±2- 5/,±3+ f^±4 + 0iiAw)') 

/^T2 = + 15/. - 40/,±i+ 45/,±2 - 24/,±3 + 5/,±4 + 0((A«;)5) 
/^T3 = + 35/, - 105/,±i + 126/,±2 - 70/,±3 + 15/,±4 + OUAwf) 

Part (b): 6 point 5th order (quintic) extrapolants: 

Ui = + 6/, -15/,±i+ 20/,±2 - 15/,±3 + 6/,±4- /,±5 + 0((At^f ) 
/,^2 = + 21/, - 70/,±i + 105/,±2 - 84/,±3 + 35/,±4 - 6f^±5 + OiiAwf) 

TABLE IV. This table gives the coefficients of our Lagrange (polynomial) extrapolants for use 
in hnite differencing. As discussed in the text, when computing d^ we always use the extrapolants in 
part (a); when computing dww we may use the extrapolants of either part (a) or (b), corresponding 
to using the molecules of part (b) or (c) of table III, respectively. 
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TABLE V. This table gives the initial data and other parameters for our test computations. 
For each model this table shows a descriptive name, the spatial grid resolution, the spatial grid 
size, the coordinate time for which we have evolved the system, the input perturbation for our 
initial data solver (cf. paper I), a qualitative description of the resulting initial data, a summary 
of the initial slice's mass distribution (broken down into black hole, scalar field, and total mass), 
and the final black hole mass at the end of the evolution. 

^ These evolutions use an alternate finite differencing scheme near the grid boundaries, as discussed 
in sections VII B 3 and IX E. 

These evolutions have a small amount of white noise added to the state vector each time the 
right-hand-side function in the evolution equations is evaluated, as discussed in appendix D. 
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Thornburg Figure 1(b) 
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Thornburg Figure 1: legend for parts (a) and (b) 
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Thornburg Figure 12: legend for parts (a) and (b) 
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> P^P + 0.02xGaussian(r'init=20,(T=5) 0.0735xGaussian(r'= 21.8 ,(T= 3.5 ) 0.976 + 0.641 = 1.617 1.183 
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